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The simplest theories often have much merit and many limitations, and in this vein, the value 
of Neutral Theory (NT) has been the subject of much debate over the past 15 years. NT was 
proposed at the turn of the century by Stephen Hubbell to explain pervasive patterns observed in 
the organization of ecosystems. Its originally tepid reception among ecologists contrasted starkly 
with the excitement it caused among physicists and mathematicians. Indeed, NT spawned several 
theoretical studies that attempted to explain empirical data and predicted trends of quantities that 
had not yet been studied. While there are a few reviews of NT oriented towards ecologists, our goal 
here is to review the quantitative results of NT and its extensions for physicists who are interested 
in learning what NT is, what its successes are and what important problems remain unresolved. 
Furthermore, we hope that this review could also be of interest to theoretical ecologists because 
many potentially interesting results are buried in the vast NT literature. We propose to make these 
more accessible by extracting them and presenting them in a logical fashion. We conclude the review 
by discussing how one might introduce realistic non-neutral elements into the current models. 
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I. INTRODUCTION 


It is interesting to contemplate an entangled hank, clothed with many plants of many kinds, with birds singing on 
the bushes, with various insects flitting about, and with worms crawling through the damp earth, and to reflect that 
these elaborately constructed forms, so different from each other in so complex a manner, have been all produced 
by laws acting around us. In this celebrated text from the Origin of Species |42j . Darwin eloquently conveys his 
amazement for the underlying simplicity of Nature: despite the striking diversity of shapes and forms, it exhibits deep 
commonalities that have emerged over wide scales of space, time and organizational complexity. For more than fifty 
years now, ecologists have collected census data for several ecosystems around the world from diverse communities 
such as tropical forests, coral reefs, plankton, etc. However, despite the contrasting biological and environmental 
conditions in these ecological communities, some macro-ecological patterns can be detected that reflect strikingly 
similar characteristics in very different communities (see BOX 1). This suggests that there are ecological mechanisms 
that are insensitive to the details of the systems and that can structure general patterns. Although the biological 
properties of individual species and their interactions retain their importance in many respects, it is likely that the 
processes that generate such macro-ecological patterns are common to a variety of ecosystems and they can therefore 
be considered to be universal. The question then is to understand how these patterns arise from just a few simple 
key features shared by all ecosystems. Contrary to inanimate matter, living organisms adapt and evolve through the 
key elements of inheritance, mutation and selection. 

This fascinating intellectual challenge fits perfectly into the way physicists approach scientihc problems and their 
style of inquiry. Statistical physics and thermodynamics have taught us an important lesson, that not all microscopic 
ingredients are equally important if a macroscopic description is all one desires. Consider for example a simple system 
like a gas. In the case of an ideal gas, the assumptions are that the molecules behave as point-like particles that do not 
interact and that only exchange energy with the walls of the container in which they are kept at a given temperature. 
Despite its vast simplifications, the theory yields amazingly accurate predictions of a multitude of phenomena, at least 
in a low-density regime and/or at not too low temperatures. Just as statistical mechanics provides a framework to 
relate the microscopic properties of individual atoms and molecules to the macroscopic or bulk properties of materials, 
ecology needs a theory to relate key biological properties at the individual scale, with macro-ecological properties at 
the community scale. Nevertheless, this step is more than a mere generalization of the standard statistical mechanics 
approach. Indeed, in contrast to inanimate matter, for which particles have a given identity with known interactions 
that are always at play, in ecosystems we deal with entities that evolve, mutate and change, and that can turn on or 
off as well as tune their interactions with partners. Thus the problem at the core of the statistical physics of ecological 
systems is to identify the key elements one needs to incorporate in models in order to reproduce the known emergent 
patterns and eventually discover new ones. 

Historically, the first models defining the dynamics of interacting ecological species were those of Lotka & Volterra, 
which describe asymmetrical interactions between predator-prey or resource-consumers systems. These models are 
based on Cause’s competitive exclusion principle |55j . which states that two species cannot occupy the same niche 
in the same environment for a long time (see BOX 2). There have been several variants and generalizations of these 
models, e.g. [32], yet all of them have several drawbacks: 1) They are mostly deterministic models and often do not 
take into account stochastic effects in the demographic dynamics [37]; 2) As the number of species in the system 
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increases, they become analytically intractable and computationally expensive; 3) They have a lot of parameters that 
are difficult to estimate from ecological data or experiments; 4) It is very difficult to draw generalizations that include 
spatial degrees of freedom; and 5) While time series of abundance are easily analyzed, it remains challenging to study 
analytically the macroecological patterns they generate and thus, their universal properties. 

A pioneering attempt to explain macro-ecological patterns as a dynamic equilibrium of basic and universal ecological 
processes - and that also implicitly introduced the concept of neutrality in ecology - was made by MacArthur and 
Wilson in the famous monograph of 1967 titled “The theory of island biogeography” [93]. In this work, the authors 
proposed that the number of species present on an island (and forming a local community) changes as the result of 
two opposing forces: on the one hand, species not yet present on the island can reach the island from the mainland 
(where there is a meta-community); and on the other hand, the species already present on the island may become 
extinct. MacArthur and Wilson’s model implies a radical departure from the then main current of thought among 
contemporary ecologists for at least three reasons: 1) Their theory stresses that demographic and environmental 
stochasticity can play a role in structuring the community as part of the classical principle of competitive exclusion; 

2) The number of coexisting species is the result of a dynamic balance between the rates of immigration and extinction; 

3) No matter which species contribute to this dynamic balance between immigration and extinction on the island, all 
the species are treated as identical. Therefore, at the level of species, they introduced a concept that is now known 
as neutrality (see BOX 2). 

Just a few years later, the American ecologist, H. Caswell, proposed a model in which the species in a commu¬ 
nity are essentially a collection of non-interacting entities and their abundance is driven solely by random migra¬ 
tion/immigration. In contrast to the mainstream vision of niche community assembly, where species persist in the 
community because they adapt to the habitat, Caswell stressed the importance of random dispersal in shaping eco¬ 
logical communities. Although the model was unable to correctly describe the empirical trends observed in a real 
ecosystem, it is important because it pictured ecosystems as an open system, within which various species have come 
together by chance, past history and random diffusion. 

Greatly inspired by the theory of island biogeography and the dispersal limitation concept (see BOX 2), in 2001 
Hubbell published an influential monograph titled “The Unified Neutral Theory of Biodiversity and Biogeography”. 
Unlike the niche theory and the approach adopted by Lotka & Volterra, the neutral theory (NT) aims to only model 
species on the same trophic level (monotrophic communities, see BOX 2), species that therefore compete with each 
other because they all feed on the same pool of limited resources. For instance, competition arises among plant species 
in a forest because all of them place demands on similar resources like carbon, light or nitrate. Other examples include 
species of corals, bees, hoverflies, butterflies, birds and so on. The NT is an ecological theory based on random drift, 
whereby organisms in the community are essentially identical in terms of their per capita probabilities of giving birth, 
dying, migrating and speciating. Thus, from an ecological point of view, the originality of Hubbell’s NT lies in the 
combination of several factors: i) it assumes competitive equivalence among interacting species; ii) it is an individual- 
based stochastic theory founded on mechanistic assumptions about the processes controlling the origin and interaction 
of biological populations at the individual level (i.e. speciation, birth, death and migration); in) it can be formulated 
as a dispersal limited sampling theory; iv) it is able to describe several macro-ecological patterns through just a few 
fundamental ecological processes, such as birth, death and migration [iHiiiniisniiTS]. Although the theory has been 
highly criticized by many ecologists as being unrealistic [351 [311 11081112311124] , it does provide very good results 
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when describing observed ecological patterns and it is simple enough to allow analytical treatment [75J 110611147H149] . 
However, such precision does not necessarily imply that communities are truly neutral and indeed, non-neutral models 
can also produce similar patterns mm]- Yet the NT does call into question approaches that are either more complex 
or equally unrealistic [10911121) . Moreover, NT is not only a useful tool to reveal universal patterns but also it is a 
framework that provides valuable information when it fails. Accordingly, these features together have made NT an 
important approach in the study of biodiversity [B) [231301 EH 1311 11117L112811152] . 

From a physicist’s perspective, NT is appealing as it represents a sort of “thermodynamic” theory of ecosystems. 
Similar to the kinetic theory of ideal gases in physics, NT is a basic theory that provides the essential ingredients to 
further explore theories that involve more complex assumptions. Indeed, NT captures the fundamental approach of 
physicists, which can be summarized by Einstein’s celebrated quote “Make everything as simple as possible, but not 
simpler”. Finally, it should be noted that the NT of biodiversity is basically the analogue of the theory of neutral 
evolution in population genetics [86] and indeed, several results obtained in population genetics can be mapped to 
the corresponding ecological case mi- 

statistical physics is contributing decisively to our understanding of biological and ecological systems by providing 
powerful theoretical tools and innovative steps [Si |S7] to understand empirical data about emerging patterns of 
biodiversity. The aim of this review is not to present a complete and exhaustive summary of all the contributions to 
this field in recent years - a goal that would be almost impossible in such an active and broad interdisciplinary field 
- but rather, we would like to introduce this exciting new field to physicists that have no background in ecology and 
yet are interested in learning about NT. Thus, we will focus on what has already been done and what issues must be 
addressed most urgently in this nascent field, that of statistical physics applied to ecological systems. A nice feature 
of this field is the availability of ecological data that can be used to falsify models and highlight their limitations. 
At the same time, we will see how the development of a quantitative theoretical framework will enable one to better 
understand the multiplicity of empirical experiments and ecological data. 

This review is organized into five main sections. Sec. [^ is an attempt to review several important results that 
have been obtained by solving neutral models at stationarity. In particular, we will present the theoretical framework 
based on Markovian assumptions to model ecological communities, where different models may be seen as the results 
of different NT ensembles. We will also show how NT, despite its simplicity, can describe patterns observed in real 
ecosystems. In Sec. we will present more recent results on dynamic quantities related to NT. In particular, we will 
discuss the continuum limit approximation of the discrete Markovian framework, paying special attention to boundary 
conditions, a subtle aspect of the time dependent solution of the NT. In Sec. |IV[ we will provide examples of how 
space plays an essential role in shaping the organization of an ecosystem. We will discuss both phenomenological, 
and spatially implicit and explicit NT models. A final subsection will be devoted to the modeling of environmental 
fragmentation and habitat loss. In Sec. |V]we will propose some emerging topics in this fledgling field, and present 
the problems currently being faced. Finally, we will close the review with a section dedicated to conclusions. 
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FIG. 1. [From [TU]]. Visual Scheme of two important macro ecological patterns (see BOX 1): RSA and SAR. The functional 
shape of the RSA depends on the spatial scale considered, while the SAR generally displays a tri-phasic behavior (see Section 


IVI. There is a growing appreciation that the various descriptors of biodiversity are intrinsically inter-related, and substantial 


efforts have been devoted to understand the links between them m- 











BOX 1: Macro-ecological patterns 


• Alpha-diversity: the number of species found in a given area, regardless of their abundance and spatial 
distributions. Sometimes it is more appropriate to assess the diversity of a local community by taking into 
account their abundances as well. For this purpose two alternative indexes have been used, the Shannon 
(H) and Simpson (D) index |132j 

• Beta-diversity: the probability that two individuals at distance r belong to the same species (i.e., they 
are conspecific). Several alternative definitions have been used in the literature, including the pair (or 
two-point) correlation function (PCF), the Sprensen similarity index (SSI) or the Jaccard similarity index 
(JSI). However, the general purpose of beta-diversity is to describe the turnover of multiple species in space. 

• Pair Correlation Function (PCF): the correlation in species’ abundance between pairs of samples as a 
function of their distance. 

• Relative Species Abundance (RSA): the probability that a species has n individuals in a given region. 
When multiplied by the total number of species in the region, this gives the number of species with n 
individuals (see section 0. This is sometimes called the Species Abundance Distribution (SAD). 

• Species (time) Turnover Distribution (STD): the probability density function that the ratio of the 
future to the current population sizes of any species has a value A in a given ecological community (see 


section III). 


Persistence or lifetime distribution (SPT): the probability density function of the time interval be¬ 


tween the emergence and local extinction of any species within a given area (see section III). 


Species Area Relationship (SAR): the function that relates the mean number of species, S', to the 
area. A, they live in. On a relatively large range of spatial scales it is well approximated by a power law, 
S(A) = cA^ where 0 < z < 1. It has been shown that the SAR in the log-log scale has three qualitatively 
different behaviours from local to continental spatial scales: approximately linear-like at very small and 
very large scales, and power-law-like for intermediate scales. This is referred to as a triphasic SAR (see 
section 


IV). 


Endemic Area Relationship (EAR): The mean number of species that are present in the area A but 


not outside it (see section IV ). 
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BOX 2: Glossary 


• Trophic level: the set of all species belonging to the same level in the food chain. Individuals of species belonging 
to higher levels feed upon those in the lower ones, while individuals belonging to the same trophic level compete for 
the same pool of resources. Neutral theory is an ecological theory for species in one specific trophic level, whereas 
Niche theory (see below) can also deal with species at different trophic levels. 

• Cause’s competitive exclusion principle: a pair of species cannot stably co-exist if they feed upon exactly the 
same resources under the same environmental conditions | 55 |. 

• Niche theory: species can stably coexist in an ecological community if their characteristics (or traits) allow them 
to specialize on one particular set of resources or environment conditions (niches) in which they are superior to 
their competitors. In other words, they occupy different niches ESI- Such niche separation is deemed to enhance 
trade-offs and facilitate co-existence, even though there is no a priori method to identify the correct niches that 
favor co-existence. The underlying rationale of the theory is Cause’s exclusion principle, and the classical model 
is Lotka-Volterra’s set of differential equations. In contrast to Niche Theory, Neutral theory claims that niche 
differences are not essential to co-existence. 

• Species-level neutral models: these models assume that all species are equivalent as they all have the same 
probability of immigration, extinction and speciation. The only state variable of these models is the number of 
species in a community and thus, it cannot predict the distribution of the population sizes across species (see RSA 
below) | 93| . 

• Individual-level neutral models: these models assume that all species are equivalent at the individual level, 
having the same birth, death, immigration and speciation rates regardless of their identity. Therefore, species 
are competitively equivalent. The state variable of these models is the population of any species in a region and 
therefore, such models can be used to understand species richness and abundance 1 75 1. 

• Symmetric models: any model whose outcomes are invariant when exchanging species identities. The family of 
symmetric models is larger than the neutral ones: for instance, effects of density dependence of individuals on the 
per capita birth and death rates can be accommodated in symmetric but not neutral models. Some authors do not 
make a distinction between symmetric and neutral models. 

• Dispersal limited process: any process that constrains offspring to disperse in the vicinity of its parents. 

• Multispecies Voter Model with Speciation (MVM): a spatially explicit neutral model in which each individual 
is located in a regular lattice and belongs to a species. Given a spatial configuration at time t, the configuration at 
time t -I- 1 is obtained as follows: an individual is chosen at random and is replaced by a copy of one of its nearest 
neighbours (chosen at random) with a probability 1 — v, ox by an individual belonging to a new species (not already 
present in the lattice) with probability v. v is the speciation parameter, although it may incorporate immigration 

_ effects. _ 
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II. NEUTRAL THEORY AT STATIONARITY 

Neutral theory deals with ecological communities within a single trophic level, i.e. communities whose species 
compete for the same pool of resources (see BOX 2). This means that neutral models will generally be tested on data 
describing species that occupy the same position in the food chain, like trees in a forest, breeding birds in a given 
region, butterflies in a landscape, plankton, etc.. Therefore, ecological food webs with predator-prey type interactions 
are not suitable to be studied with standard neutral models. 

As explained in the introduction, ecologists have been studying an array of biodiversity descriptors over the last 
sixty years (see BOX 1), including relative species abundance distributions (RSA), species-area relationships (SAR) 
and spatial pair correlation function (PCF). For instance, the RSA represents one of the most commonly used static 
measures to summarize information on ecosystem diversity. The analysis of this pattern reveals that the RSA distri¬ 
butions in tropical forests share similar shapes, regardless of the type of ecosystem, geographical location or the details 
of species interactions (see Fig. [^. Therefore, the functional form of the RSA (see seminal papers by Fisher and 
Preston |118j for a theoretical explanation of its origins) has been one of the great problems studied by ecologists. 
Indeed, a great deal of attention has been devoted to the precise functional forms of these patterns. 

It is instructive to start deriving some of these within a very simple but extreme neutral model, which assumes that 
species are independent and randomly distributed in space. This null model tells us what we should expect when the 
observed macroecological patterns are only driven by randomness, with no underlying ecological mechanisms. If the 
density of individuals in a very large region is p, then the probability that a species has n individuals within an area 
a is well approximated by a Poisson distribution with a mean pa. As defined in Box 1, this is the RSA for the area a. 
However, the empirical data are not well described by a Poisson distribution and as we shall see later on, better fits 
are usually given by log-series, gamma or log-normal distributions. 

The SAR curve can also be calculated within a slightly more accurate model, which still assumes that species are 
independent and randomly situated in space. Let us now suppose that a region with area Aq contains Stot species in 
total (of-diversity - see Box 1), and that the species i has rii individuals in Aq. If we consider a smaller area A within 
the region, then the probability that an individual will not be found in such area is 1 — A/Aq, while the probability 
that the whole species i is not present therein is (1 — A/Aq)”' = 1—pi. If we now consider the random variable Ii{A), 
which is 1 when the species i is found within the area A and 0 if not, then (/^(A)) = = 1 — (1 — A/Ao)"S because 

Ii{A) is a Bernoulli random variable with an expectation value Pi- Therefore, the mean number of species in the area 
A (i.e. the SAR) is simply S'(A) = which is 


Stot 

S{A) = Stt,t-Y.^l-A/Aor. ( 1 ) 

i=l 

Although this model was originally studied by Coleman [39], we now know that it significantly overestimates species 
diversity at almost all spatial scales m- 

Beta-diversity (see Box 1) can be estimated under the assumptions we have mentioned. Now, regardless of the spatial 
distance between two individuals, the probability that two of them belong to the same species i is ni{ni — l)/[N{N—l)], 
where N is the total number of individuals in the community, i.e. N = Therefore, the probability to find any 

pair of con-specific individuals is { 1 / 5 *°*) X]i=i ~ l)/[-^(-^~ !)]• This means that the random placement model 
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with independent species predicts that beta-diversity should not depend on the distance between two individuals. 
Again, we now have clear evidence that the probability that two individuals at distance r belong to the same species 
is a decaying function of r |104j . 

The failure of the random placement model to capture the RSA, SAR and beta-diversity is a clear indication that 
ecological patterns are driven by non-trivial mechanisms that need to be appropriately identified. Thus, we shall 
assess to what extent the NT at stationarity can provide predictions in agreement with empirical data. 
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FIG. 2. [From | 149 |1. Tree relative species abundance data from the BCI, Yasuni, Pasoh, Lambir, Korup, and Sinharaja plots, 
for trees that are 10 cm in stem diameter at breast height. The frequency distributions are plotted using Preston’s binning 
method as described in |148| and the bars are the observed number of species binned into log(2) abundance categories. The 
first histogram bar represents the second bar -|- the third bar -|- (c^s) -|- the fourth bar -|- ((^ 5 ) -|- 

{4>6) + {4>r) + 8 -iid so on. Here {ij>n) is the number of species with an abundance n. As examples, we show the fits of a 
density-dependent symmetric model (red line) and dispersal limitation model (blue circles), which will be studied in greater 
detail in sec. [n|and[rv] respectively. 


There are two related, but distinct analytical frameworks that have been used to mathematically formulate the NT of 
biodiversity at stationarity for both local and met a-communities. From an ecological point of view, a local community 
is defined as a group of potentially interacting species sharing the same environment and resources. Mathematically, 
when modeling a local community the total community population abundance remains fixed. Alternatively, a meta¬ 
community can be considered a set of interacting communities that are linked by dispersal and migration phenomena. 
In this case, it is the average total abundance of the whole meta-community that is held constant. From the physical 
point of view, these roughly correspond to the micro/canonical (fixed total abundance) and the grand canonical (fixed 
average total abundance) ensembles, respectively. The micro/canonical ensemble or so-called zero sum dynamics 
when death and birth events always occur as a pair originates from the sampling frameworks in population genetics 
pioneered by Warren Evens and Ronald Fisher [^. It should be noted that even though a fixed size sample is one 
way to analyze available data, for the majority of cases (apart from very small size samples), the grand canonical 
ensemble approach is that used routinely in statistical physics and it provides a very precise yet largely simplified 
description of the system. The ultimate reason for this lies in the surprising accuracy of the asymptotic expansion 
of the gamma function (the mathematical framework heavily uses combinatorials and factorials etc.). The Stirling 
approximation can be used for very large values of the gamma function, nevertheless, it is quite accurate even for 
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values of the arguments of the order of 20. The advantages of the master equation (see below) and the grand canonical 
ensemble approach stem from their computational simplicity, which make the results more intuitively transparent. 

We now introduce the mathematical tools of stochastic processes that will be used extensively in the rest of the 
article. 


A. Markovian modeling of neutral ecological communities 


1. The Master Equation of birth and death 


Let c be a configuration of an ecosystem that could be as detailed as the characteristics of all individuals in the 
ecosystem, including their spatial locations or as minimal as the abundance of a specified subset of species. Let 
P(c, <|co,to) be the probability that a configuration c is seen at time t, given that the configuration at time to was 
Co (referred to as P{c,t) for simplicity). For our applications, the configurations c are typically species abundance 
(denoted by n). 

Assuming that the stochastic dynamics are Markovian, the time evolution of P{c, t) is given by the Master Equation 

(ME) [MiiHg 

('r[c|c']P(c',t)-T[c'|c]P(c,t)Y (2) 

^/ \ / 

where T[c|c'] is the transition rate from the configuration P to configuration c. Under suitable and very plausible 
conditions [53], P{c,t) approaches a stationary solution at long times, Ps(c), which satisfies the following equation 

^ (t[c\c']Ps{P) - T[P\c]P,{c)] = 0 . ( 3 ) 

c' ^ ' 

This equation is typically intractable with analytical tools, because it involves a sum of all the configurations. If each 
term in the summation is zero, i.e. 

r[c|c']p«(c')-r[c'|c]p,(c) = o, (4) 

detailed balance is said to hold. A necessary and sufficient condition for the validity of detailed balance is that for 
all possible cycles in the configuration space, the probability of walking through it in one direction is equal to the 
probability of walking through it in the opposite direction. Given a cycle {ci, C 2 ,..., c„_i, ci}, detailed balance holds 
if and only if, for every such cycle, 

T[ci|c 2 ]T[c 2 |c 3 ] . . .T[c„_ 2 |c„_i]r[c„_i|ci] = P[ci | C„_ i]r[c„_ i |c„_ 2 ] . • . T[c 3 | C 2 ]r[c 2 | Cl] . (5) 


This condition evidently corresponds to a time-reversible condition. 

Now let us apply these mathematical tools to the study of community dynamics that is driven by random demo- 
graphical drift. Consider a well-mixed local community. This is equivalent to saying that the distribution of species 
in space is not relevant, which should hold for an ecosystem with a linear size smaller than, or of the same order as 
the seed dispersal range. In this case one can use c = (rii,..., ng) = n where is the population of the f-th species 
and the system contains S species. We can rewrite eq. [^in the following way: 


i9P(n, t) 


T[n\n']P{n\t) - ^ P[n'|n]P(n, t). 

n'^n n'^n 


dt 


(6) 
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In this case, T[n'|n] can take into account birth and death, as well as immigration from a meta-community. The 
simplest hypothesis is that T[n'|n] is the result of S elementary birth and death processes that occur independently 
for each of the S species i.e. 

s 

^[n'|n] = X! n {Sn'^,nk+iTk{n + l\n) + - l|n)^ , (7) 

k—1 i^k 


where 


Tk{n + l|n) = b{n, k) 

is the birth rate and 


( 8 ) 


Tkin - l|n) = d(n, fc). 


(9) 


the death rate. This particular choice corresponds to a sort of mean-field (ME) approach [71 1103111121114011148111491 

dMl. 

Our many-body ecological system can also be formulated in a language more familiar to statistical physicists, where 
we consider the distribution of balls into boxes. The “boxes” are the species and the “balls” are the individuals. 
Birth/death processes correspond to adding/removing a ball to/from one of the boxes using a rule as dictated by Eqs. 
([^ and (H). Eq. (§ with the choice 0 can be simplified if one assumes that the initial condition is factorized as 
P{n,t = 0) = nfe=i Pk{nk-,t = 0). In this case the solution is again factorized as P{n,t) — 11^=1 Pk{nk,t) where 
each Pk satisfies the following ME in one degree of freedom Uk 

dPk{nk:t) ^ Pk{rik + l,t)d{nk + 1, A:) - {b{nk,k) + d{nk,k)) Pk{nk,t) . (10) 


The stationary solution of eq.(lO) is easily seen to satisfy detailed balance [oil 183] . 

First, we note that because of the neutrality hypothesis, species are assumed to be demographically identical and 
therefore, we can drop the k dependence factor from the equations @-(|9|. In other words, we can concentrate on the 
probability P{n) that a given species (box) has an abundance n (balls). In this case, species do not interact and thus, 
in our calculation we can follow a particular species (the boxes are taken to be independent). Following equation Q, 
it is easy to see that the solution of the birth-death ME 0 that satisfies the detailed balance condition and that thus 
corresponds to equilibrium, is jSH IHH] 


=^on 


d{z) 


( 11 ) 


where Pg can be calculated by the normalization condition Psi'^k) = 1, and assuming that all the rates are positive. 
In particular, 6(0) > 0 and d{0) = 0. 

When there are S boxes, all satisfying the same birth-death rules, the general and unique equilibrium solution is 

s 

Psini,n2,...,ns) = Y[Ps{nk) ( 12 ) 

k=l 

Depending on the functional form of b{n) and d{n), one can readily work out the desired Pg{n). We will start with 
some cases that are familiar to physicists |I51j , and then move onto more ecologically meaningful cases |148[ 1149] . 
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2. Physics Ensembles 


The Random walk and Bose-Einstein Distribution. If one chooses b{n) = bo and d{n) = di for n > 1, and 
d{n) = 0 otherwise (or alternatively, b{n) = (n + l)6o ^^nd d{n) = din), one obtains a pure exponential distribution 
P{n) = r"(l — r) where r = bo/di. Note that in both these cases, b{n — l)/d{n) is the same. Substituting this in 
equation ( [l^ gives us the well known Bose-Einstein distribution for non-degenerate energy levels, i.e. 

P(ni,n 2 ,...,ns) = r^(l-r)‘^, (13) 

where N = Here, r corresponds to in the grand canonical ensemble. 

The Fermi-Dirac Distribution. If d{n) = din and b{n) = 0 for any n other than 0, and equal to bg otherwise, 
then P{n) = r"/(l -I- r) for n = 0 or n = 1 and P{n) = 0 for other values of n, and accordingly the Fermi-Dirac 
distribution is achieved 


P(ni, 712,..., ns) = "'*=(!-|-r) ^ for = 0 or = 1 Vz (14) 

P(ni, 772,..., Tis) = 0 otherwise. (15) 


Boltzmann counting If d{n) = din and = bo for all n, then one obtains a Poisson distribution P{n) = e '^r'^/nl, 
and this leads to 


J.N 

P(77i, 772,..., 77s) OC —- -. (16) 

\Ak=ink'- 

This is the familiar grand-canonical ensemble Boltzmann counting in physics, where r plays the role of fugacity. It is 
noteworthy that, unlike the conventional classical treatment iza, where an additional factor of N\ is obtained, here 
one gets the correct Boltzmann counting and thereby avoids the well-known Gibbs paradox [74]. Thus, if one were 
to ascribe energy values to each of the boxes and enforce a fixed average total energy, the standard Boltzmann result 
would be obtained whereby the probability of occupancy of an energy level e is proportional to e~^^, where (3 is 
proportional to the inverse of the temperature. 


3. Ecological Ensembles 

Density independent dynamics. We now consider the dynamic rules of birth, death and speciation that govern 
the population of an individual species. The most simple ecologically meaningful case is to consider d{n) = dn and 
b{n) = bn for ti > 0, and 6(0) > 0, d{0) = 6(— 1) = 0. We now define r = b/d. Moreover, in order to ensure that the 
community will not become extinct at longer times, speciation may be introduced by ascribing a non-zero probability 
of the appearance of an individual from a new species, i.e. 6(0) = ho = v. In this case the probability for a species of 
having n individuals at stationarity is: 

P(r7) = i>[l — £'ln(l — r)]“^r"/77, (17) 

where n = v/h and 77 > 0 and P(0) = 1/(1 — — r)). 

The RSA {(j){n)) is the average number of species with a population n, and this is simply |148j 


{(j){n)) = SP{n) = dr'^jn. 


(18) 
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This is the celebrated Fisher log-series distribution, i.e. the distribution Fisher proposed as that describing the 
empirical RSA in real ecosystems [52] . The parameter 9 = Sv/[1 — i/ln(l — r)] is known as the Fisher number or 
biodiversity parameter. 

In 1948, another great ecologist of the inductive school, F.W. Preston, published a paper m challenging Fisher’s 
point of view. He showed that the Log-series is not a good description for the data from a large sample of birds. In 
fact, he observed an internal mode in the RSA that was absent in a Log-series distribution. In particular, Preston 
introduced a way to plot the experimental RSA data by octave abundance classes (i.e. [2^, 2^+^], for fc = 1, 2,3,...), 
showing that a good fit of the data was represented by a Log-Normal distribution. Indeed, there are several examples 
of RSA data that display this internal mode feature (see Fig. [^. The intuition of Preston was that the shape of 
the RSA must depend on the sampling intensity or size of the community. Conventionally, when studying ecological 
communities, ecologists separate them into two distinct classes: small local communities (e.g., on a island) and meta¬ 
communities of much larger communities or those composed of several smaller local communities (see Fig. [^. The 
neutral modeling schemes for these two cases - that we will denote by the sub-script L and M respectively - are 
not the same as the ecological processes involved differ. In fact, the immigration rate (to) in a local community 
is a crucial parameter as the community is mainly structured by dispersal limited mechanisms, and the speciation 
rate {v) can be neglected. On the other hand, in a meta-community immigration does not occur (species colonize 
within the community) and the community is shaped by birth-death-competition processes, although speciation v also 
plays an important role. A meta-community can also be thought of as consisting of many small semi-isolated local 
communities, each of which receives immigrants from other local communities. When considering meta-community 
dynamics, the natural choice is to put a soft constraint on Jm, he. the total number of individuals is free to fluctuate 
around the average population {Jm)- Indeed, one finds that at the largest Jm limit, the results obtained with a 
hard constraint on Jm are equivalent to those with a soft constraint. On the other hand, when considering a local 
community, it is safer to place a hard constraint on the total population J^. In both these cases, S represents the total 
number of species that may potentially be present in the community, while the average number of species observed 
in the community is denoted by {S). 

There are several ecological meaningful mechanisms that can generate a bell shaped Preston-like RSA. The first of 
these involves density dependent effects on birth and death rates. The second involves considering a Fisher log-series 
as the RSA of a meta-community acting as a source of immigrants to a local community embedded within it. The 
dynamics of the local community are governed by births, deaths and immigration, whereas the meta-community 
is characterized by births, deaths and speciation. This leads to a local community RSA with an internal mode 
[102111401114811149] . A third way is to incrementally aggregate several local communities (see Appendix A) 

Local Dynamics with Density Dependent Birth Rates One major puzzle in community ecology is the 
coexistence of a large number of tree species on a local scale in tropical forests. This phenomenon is often explained 
by invoking density- and frequency dependent mechanisms. Processes that hold the abundance of a common species 
in check inevitably lead to rare-species advantages, given that the space or resources freed up by density-dependent 
death can be exploited by less-common species. Therefore, inter-species frequency dependency is the community- 
level consequence of intra-species density dependence, and thus, they are two different manifestations of the same 
phenomenon [147] . 
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We begin by noting that the mean number of species with n individuals, (</'«), is not determined by the absolute rates 
of birth or death but rather, by their ratio, = d(i+i\) • This follows from the observation that (0„) is proportional 
to {ri,kr 2 ,k ■ ■ ■‘rn-i,k)k, where the average {■ ■■)k is obtained from all the species. This simple formulation |147) is 
sufficiently general to represent the communities of either symmetric species (in which all the species have the same 
demographic birth and death rates) or the case of asymmetric or distinct species. The more general asymmetric 
situation captures niche differences and/or differing immigration fiuxes that might arise from the different relative 
abundances of distinct species in the surrounding meta-community (see BOX 2). 

Two of the most prominent hypotheses to explain frequency and density dependence are the Janzen-Connell 1111179] 
and the Chesson-Warner hypotheses |36j . These mechanisms generally predict the reproductive advantage of a rare 
species due to ecological factors and they can be readily captured in a common mathematical framework that will be 
presented below. 

The Janzen-Connell hypothesis postulates that host-specific pathogens or predators act in the vicinity of the ma¬ 
ternal parent. Thus, seeds that disperse further away from the mother are more likely to escape mortality. This 
spatially structured mortality effect suppresses the uncontrolled population growth of locally abundant species rela¬ 
tive to uncommon species, thereby producing reproductive advantage to a rare species. The Chesson-Warner storage 
hypothesis explores the consequences of a variable external environment and it relies on three empirically validated 
observations: species respond in a species-specific manner to the fluctuating environment; there is a covariance be¬ 
tween the environment, and intra- and inter-species competition in function of the abundances of the species; life 
history stages buffer the growth of population against unfavorable conditions. Such conditions prevail when species 
have similar per capita rates of mortality but they reproduce asynchronously and there are overlapping generations. 

We now introduce a modified symmetric theory that captures density- and frequency dependence (rare species 
advantage or common species disadvantage), and in which , should be a decreasing function of abundance in order 
to incorporate density dependence. The equations of density dependence in the per capita birth and death rates for 
an arbitrary species of abundance n are: = b ■ [l -|- ^ -|- o (;^)] and = d ■ [l -|- ^ -|- o (;^)], for n > 0 as 

the leading term of a power series in (1/n), = b ■ ~ where bi and di are 

constants. This expansion captures the essence of density-dependence by ensuring that the per-capita rates decrease 
and approach a constant value for a large n, given that the higher order terms are negligible. As noted earlier, the 
quantity that controls the RSA distribution is the ratio bn/dn+i- Thus, the birth and death rates, and dm can be 
defined up to /(n -|- 1) and f{n) respectively, where / is any arbitrary well-behaved function. 

Strikingly, any relative abundance data can be considered as arising from effective density dependent processes in 
which the birth and death rates are given by the above expressions. Thus, one would expect that the per capita birth 
rate or fecundity drops as the abundance increases, whereas mortality ought to increase with abundance. Indeed, 
the per capita death rate can be arranged to be an increasing function of n, as observed in nature, by choosing an 
appropriate function / and appropriately adjusting the birth rate so that the ratio bnjdn+i remains the same. This 
then ensures that the RSA does not vary. 

The mathematical formulation of density dependence may seem unusual to ecologists familiar with the logistic or 
Lotka-Volterra systems of equations, wherein density dependence is typically described as a polynomial expansion 
of powers of n truncated at the quadratic level. Such an expansion is valid when the characteristic scale of n is 
determined by a fixed carrying capacity. Conversely, here the range of n is from 1 to an arbitrarily large value and 





17 


not to some carrying capacity. Therefore, an expansion in terms of powers of (1/n.) is more appropriate. For this 
symmetric model |147j . and bearing in mind that {4>n) = SPq rir=/ ; one readily arrives at the following relative 

species-abundance relationship: 


{K) = 0 


n + c' 


(19) 


where x = b/d and for parsimony, we make the simple assumption that bi = di = c and the higher order coefficients, 
^ 2 , ^ 2 , 63 , da, ..., are all 0. The biodiversity parameter, 0, is the normalization constant that ensures the total 
number of species in the community is S, and it is given hy 6 = S^^/F{1 + c,2 + c,x), where F{l-\-c,2 + c,x) is the 
standard hypergeometric function. The parameter c measures the strength of the symmetric density- and frequency 
dependence in the community, and it controls the shape of the RSA distribution. This simple model |147j does a 
good job of matching the patterns of abundance distribution observed in the tropical forest communities throughout 
the world (see Fig. [^. Note that when c —)■ 0 (the case of no density dependence), one obtains the Fisher log-series. 
In this case 9 captures the effects of speciation. 


Local community with immigrants from a meta-community. Various analytical solutions for the Markovian 
(ME) approach have been suggested in the literature. Based on pioneering work by McKane, Alonso and Sole |102) . 
an analytical solution was proposed for the ME when the size of both the local community and the meta-community 
was fixed [TTMiTiniim] . Here, we will refer to the work of Vallade and Houchmandzadeh mni, in whose model the 
birth and death rates in the meta-community are given by 


T{n + l|n) = bM{n) = 
T{n — l|n) = dM{n) = 


Jm — n 


n 


-( 1 -^) 


Jm Jm — 1 

n[JM — n + {n — l)v\ 

Jm{Jm — 1 ) 


( 20 ) 

( 21 ) 


Eq. (20) describes the event that a randomly chosen individual of the species under consideration (probability n/ Jm) 
colonizes (probability 1 — ^) a site occupied by a different species (probability {J — M — n)/J m- Alternatively, Eq. 


( 21 ) gives the probability that this species is removed (probability n/ Jm) and replaced by another species from within 


(with probability {l — v){JM — n)/{JM — ^)) or outside (migration event with probability v) the system. Eqs. (20l-(21) 
also correspond to the transition rates for the mean field Voter Model (see MVM Box 2 and Section IV). 

Let {4>n{t))M designate the average number of species of abundance n in the meta-community. Then, in the 
stationary limit (t —> 00 ), the species that contribute to {4>n{t))M are those which entered the system by mutation at 
time t — T and that have reached size n at time t |135L I140| . As speciation is a Poissonian event (of rate v) and due 
to neutrality, all species have the same probability p{T)dT = vdr of appearing between the time interval [t, r -I- dr]. 
Thus, the time evolution of {(j)n{t)) is then simply given by {(j)n{t)) = P{n,t — T)p{T)dT = u f*P(n,T)dT, where 
P{n,t) is the solution of the ME (|^ with transition rates given by equations (20 1 - © and with the initial condition 
P{n,0) = 6n,i (<5 represents the Kronecker delta). Using Laplace transforms, the stationary RSA for the meta¬ 
community of size Jm is m- 


. ^ _ 9T{Jm + l)r(JM 9 — k) 

~ fcr(jM + i-fc)F(jM + 0)’ 


( 22 ) 


where 9 = (Jm — l)j^/(l — J^) ■ 
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Therefore the density of species of relative abundance w = k/JM in the meta community can be written as: 


gM(w) 


6»(1 - w)®-! 
cu 


(23) 


where gmi^) = JM{4'k) m and with w and 0 kept hxed (implying that v tends to zero and k tends to 

infinity). 

Let us now consider a local community (of size Jl) in contact with the aforementioned meta-community. As 
Jm ^ Jl, then mutations within the local community can be neglected (v —>■ 0) and the equilibrium distribution 
in the met a-community is negligibly modified by migration from or towards the community. Nevertheless, there is 
migration from the meta-community to the local community and thus, a migration rate m must be included in the 
transition rates for the local community dynamics: 


bL{n) 

dL(n) 


Jb-n 

Jl 


— (l-m) + 


n 

Jl 


Jl -n 
Jl-1 


(1 — m) + m(l — Lu) 


(24) 

(25) 


where u = k/JM is the relative abundance of the species in the meta-community and k is the abundance of the 
tracked species within the met a-community. 

The stationary solution Pnioj) of the corresponding ME can be calculated exactly using Eq. considering m as 

the immigration rate |140j : 


n O f igu})n[f^(l - w)] 
tnfWj — I I- - r • 

\n J ( m ) ji . 


(26) 


where (■^) is the binomial coefficient, (a)„ = r(a-l-n)/r(a) is the Pochammer coefficient, and g = (Jl — l)m/(l —m). 

Finally, the average number of species with abundance n in the local community, {4>n)L, is given by the following 
expression: 


Jm 

{4>n)L =''^Pn{k/JM){(t>k)M, (27) 

k=l 

which displays an internal mode for appropriate values of the model parameters. A generalization of this result to 
a system of two communities of arbitrary yet fixed sizes that are subject to both speciation and migration has 
also been carried out. 

Alonso and McKane [7] suggested that the log-series solution for the species abundance distribution in the meta¬ 
community is applicable only for species-rich communities, and that it does not adequately describe species-poor 
meta-communities. An application of the ME approach to the asymmetric species case was considered [S] and it 
was further demonstrated that the species abundance distribution has exactly the same sampling formula for both 
zero-sum and non-zero-sum models within the neutral approximation [481 I62j . The simplest mode of speciation 
- a point mutation - has been the one most commonly used to derive the species abundance distribution of the 
met a-community. A Markovian approach incorporating various modes of speciation such as random fission was also 
presented by Haegeman and Etienne [531 El] ■ 













19 



C) Coalescent Approach 


D) Merging several Local Communities 



FIG. 3. Four distinct neutral models of community ecology. In all these models, rii represents the abundance of the species i in 
the community, and Jl and Jm the total abundance in the local and meta-communities, respectively. (A) Hubbell’s zero-sum 
neutral model In the local community, each death is immediately followed by a birth or an immigration event. Speciation 
(or, equivalently, immigration) enables diversity to be maintained in the meta-community. (B) Local community with immigrants 
from a meta-community m The local community now interacts with the meta-community through a migration process (m). 
(C) Coalescent-type approach m where community members are traced back to the ancestors that once immigrated into 
the community (D) Joint RSA of many local communities \149'l - The whole meta-community RSA distribution is built by 
considering the joint RSA distributions of multiple local communities 


B. Coalescent-type approach of Neutral Theory 

An exciting approach pioneered by Etienne and OUT, as well as by Alonso, considers a coalescent-type phenomenon 
where community members are traced back to their ancestors that once immigrated into the community [461149] . A 
succinct comparison of the two different approaches adopted is presented in [an Hz]. 

Within this framework, the local community consists of a fixed number of individuals J that undergo zero-sum 
dynamics. At each timestep a randomly chosen individual dies and is replaced by a new individual that is an offspring 
of a randomly chosen individual in the community (of size J — 1), or by one of I potential immigrants from the 
meta-community. The immigration rate is thus to = J/(/-|-J—1), where the parameter I is merely a proxy for the 
immigration rate. The meta-community is governed by the same dynamics, with speciation replacing immigration. 

Two observations are crucial to derive the sampling formula for the species abundance distribution. Firstly, because 
there is no speciation in the local community, each individual is either an immigrant or a descendent of an immigrant. 
Thus, the information pertaining to the species abundance distribution is specified by considering “the ancestral 
tree” (tracing back each individual to its immigrant ancestor - which is somewhat similar to the phylogenetic tree 
construction in genetics) and the species composition of the set of all the ancestors. Secondly, this set of all ancestors 
can be considered as a random sample from the meta-comunity and its species abundance distribution is provided by 
the well-known multivariate Ewens distribution, corresponding to the Fisher’s log-series in the limit of large sample 
size, that describes neutral population genetics models with speciation and no immigration. As we do not know the 
ancestor of each individual, the final formula for the species-abundance composition of the data comes from summing 
the probabilities of all possible combinations of ancestors that give rise to the observed data. 










20 


In a tour de force calculation, the resulting sampling formula is found to be m 

j 


P[D\9,m,J] = 




(28) 


rii^i nj=i ,4=5 

where P[D\9,m, J] is the probability of observing the empirical data D consisting of S species with abundances rii, 
i S [1) >5'] given the value of the fundamental biodiversity number 9 of the meta-community, the immigration rate m 
and the local community size J. <I>j is the number of species with abundance j, 

Gin- n \ GI n ■ II 

(29) 


KiD,A)= n 

{al,...,as| ai=A *=1 


s{ni, ai)s{ai, 1 ) 


s(ni, 1 ) 


where s is the unsigned Stirling number of the first kind. Finally, (x)y is the Pochhammer symbol 11^=1 -I- i — 1). 

These expressions are simplihed versions of the original results presented by the authors and as expected, they 
provide a comparable fit to the ME approach for the data from the tropical forests. It is heartening that equivalent 
results can be derived using very different approaches. 


III. DYNAMICS IN NEUTRAL THEORY 


So far we have focused on the implications of neutral theory when models describing neutral patterns reach a steady 
state. The stationary condition allows one to take advantage of a variety of different mathematical techniques to 
obtain analytical expressions of ecological patterns. However, stationarity is not always a good approximation, either 
because the ecosystem is still in a state of flux, or because the assumption may hide different and important processes 
that lead to the same final steady state. Furthermore, one can calculate time dependent correlation functions [SS] 


(see, for example, (46) below). It is therefore essential to understand the temporal behavior of ecosystems in order 
to discriminate between ecological processes that would otherwise be indistinguishable at stationarity. Statistical 
comparisons between time-dependent patterns are usually more difficult than those between stationary patterns 
because they require more data and long empirical time series that are rarely available. In addition, although time- 
dependent solutions facilitate stronger tests when confronted with data, they are more difficult to obtain and this is 
the reason why only a few studies have investigated the temporal behavior of neutral models. 

An important method to study the time dependence is van Kampen’s system size expansion. Assuming that the 
total population of the system, A, is very large and the initial population is of the order of A, one expects that, at 
early times at least, the probability distribution of the system, P(n,t), would peak sharply around the macroscopic 
value n = with a width of order ~ '/N. The function is chosen to describe the evolution of the peak. These 

considerations lead to the Ansatz n = N(f>{t) + \/A^, where ^ is a new random variable. When expressing the ME in 
terms of ^ and expanding it in a/A, at order -s/A one recovers the deterministic evolution of the macroscopic state of 
the system, and at order A*^ one obtains a linear Fokker-Planck equation whose coefficients depend on time through 
the function Therefore, fluctuations around the peak are Gaussian in a first approximation (details of the method 
can be found in Ref. [53]). This expansion has some attractive features but it is usually only a good approximation 
of the temporal evolution of the original ME at short or intermediate temporal scales, since the Gaussian behavior is 
usually lost over longer periods. The time dependence of mean-held neutral models was investigated by means of van 
Kampen’s system size expansion in [102] and [103] and a good agreement with simulations was found at early times. 
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Some neutral models can also be formulated in terms of discrete-time Markov processes that have the form = 

QijPjW^ where Pi{t) is the probability that a species has i individuals at time t, and Qij is the probability that 
a species changes its abundance from j to i in one time step (i.e. the transition matrix). If the community has N 
individuals in total and either one birth or one death is allowed in one time step, then Q is a tridiagonal matrix with 
dimension TV -|- 1. Solving these models in time is basically equivalent to finding the eigenvalues along with the left 
and right eigenvectors of the matrix Q, which is usually a non-trivial task. Using a theorem derived in population 
genetics, Chisholm [37] was able to obtain the eigen decomposition of the transition matrix of a neutral model that 
describes the behavior of the local community in the original Hubbell model m- Also, he was able to show that the 
first two eigenvectors are sufficient to provide a good approximation of the full time-dependent solution. 

Other authors have focused on specihc temporal patterns that can sometimes be calculated more easily. For 
instance, a measure that has been used to quantify the diversity of a community is the Simpson index [132] . D, which 
is dehned as the probability that two individuals drawn randomly from a well-mixed community belong to the same 
species. If the individuals are drawn with replacement, then one obtains D = where rii is the number of 

individuals of species i (we have already encountered this index in Section II when we calculated the beta-diversity 
in the random placement model). This index is also a good indicator of the relative importance of spatial effects 
in a community. The dynamic evolution of D is usually simple enough to allow analytical calculations [142] . Other 
interesting dynamic properties of neutral models have been studied, where species’ extinction and monodominance 
were investigated by looking into the first-passage properties and fixation times of the process [nidi]. Empirical 
power-law relationships between the temporal mean and variance of population fluctuations (the so-called Taylor’s 
power law) are also in good agreement with neutral predictions [84] . 

Although it is not mandatory from the neutral assumption, models often assume that demographic stochasticity 
is the main source of fluctuations in stochastic neutral dynamics [nl 1122] . Demographic randomness originates from 
the intrinsic stochastic nature of birth and death events within a discrete population of individuals. However, other 
sources may be important, such as environmental stochasticity that, by contrast, encompasses effects of abiotic and 
biotic environmental variables on all individuals. There is theoretical and empirical evidence that the dynamics of the 
more common species may be driven by environmental stochasticity [89] . Consequently, neutral models only based 
on demographic stochasticity will overestimate the expected times to extinction for abundant species [123] . whose 
temporal fluctuations will also be underestimated [107] . Incorporating an environmental source of randomness as well 
as more realistic forms of speciation have made some dynamic aspects of NT more realistic @| . 

Here we will focus on the dynamic aspects of the species abundance distribution (SAD) [T^j which, under appropriate 
assumptions, can be studied in detail [911196] and whose predictions can be benchmarked against empirical data. We 
will also review recent progress in modeling dynamic patterns, including the species turnover distribution m as well 
as species’ persistence-times (or lifetimes) distributions [T^ 11131111511134] (see BOX 1). 


A. The continuum limit of the Master Equation 

The microscopic description of a system starts by correctly identifying the variables that define all the possible 
configurations of the system. Having decided how to describe the states of the system, the next step is to consider the 
transition rates among different states. The configurations of an ecological community can be described by different 














22 


variables according to different levels of coarsening of the spatio-temporal scales. However, for the sake of simplicity, 
we will focus on ecosystems comprising S species and N individuals in total, in which the configurations are specified 
fully by the variable n = (ni,... ,ns), where rii is the population of the i-th species. Therefore, for the time being, 
we will ignore any spatial effects, which will be considered in section [Tvl 

If the ecological community at time <o was in the configuration ng, then the probability that the system will be 
in configuration n at time t > \b P(n, t|no,to) (-P(n, t) for brevity). Assuming a Markovian stochastic dynamics, 
the evolution of P{n,t) is governed by the ME ([^ with transition rates T[n|n'] from the states n' to n. As seen in 
section II, the specification of a model is tantamount to defining the transition rates. This equation can be recast in 
another form that allows a simple expansion. Interactions among species change the variable n in steps of size r^, 
where ^ = 1,2,..., L and L is the total number of distinct possible changes or reactions. For instance, if a species can 
either decrease or increase its population by fc in a given time step, then L = 2. Species that change by one or that 
remain constant in a given time interval are instead described by L = 3, and when the total population is conserved, 

L 

we have ^ = 0. Of all the possible reactions, however, only a few are usually significant, ecologically relevant or 

i.=i 

meaningful. Moreover, unlike chemical reactions these meaningful reactions can be considered irreversible. Therefore, 
instead of summing over states n', as in eq.(|^, we can sum over L different reactions. If t^[n] is the transition rate of 
the £-th reaction which involves the jumps when species’ populations are n, then we can recast eq. (§ in the form 


dP(n, t) 
Wt 


= ^ {ti[n - V(\P{n - rg,t) - U[n]P{n,t)} . 


(30) 


e=i 


Assuming that ti[n] is only a function of x = n/N (this corresponds to eqs.(20), (21) with N = Jm and (24), (25) 
with N = Jl and A^ ^ 1 in both cases), we can write 


9P(x,t) 

Wt 


= H 


Ni 


p{y^- -t^[x]P(x,t)} 


(31) 


where n has been replaced by Nx. This form of the ME suggests that \ri\/N <C 1 when A^ ^ 1 and hence, a Taylor 
expansion around x of ti and P is legitimate, at least far from the boundaries. This is known as the Kramers-Moyal 
expansion [54) . Truncating to second order in N, we obtain a Fokker-Planck equation of the form 


aP(x,T) 

5^ 


9 r /I / \ n/ M 1 
= - E ^ [^.(x)P(x, r)] + — ^ 


i=l 


2N ^ dxidxj 


[Pij-(x)P(x,T)]-hO(Af ) 


where the time has been rescaled (r = t/N), the functions Ai(x) have been defined as follows 


(32) 


Ai(x) = [x] 


(33) 


and the entries of the matrix P(x) are 


Bij{x) = [x] . 


e=i 


(34) 


The matrix P(x) is a A x S' symmetric matrix that is positive semi-definite |54j . For one-step jumps, the entries 


of the matrices 0 vi can only be 1, -1 or 0. Unlike the ME (30), the Fokker-Planck equation (32) governs the 
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evolution of continuous stochastic variables that represent the relative species’ populations under the so-called diffusion 
approximation. As the total population extends to infinity, we recover the deterministic evolution given by the Liouville 
equation 




dr ^ dxi 

i—1 


(35) 


whose solution has the simple form P(x, r) = i5(x — x(r)), where 5 is the Dirac delta function. This form of P implies 
that x(r), i.e. the relative species’ populations, are the solutions of the deterministic system Xiir) = Ai(x), where 

i = 1,2,..., S [SI]- 

So far we have assumed that the system has a fixed and finite total population, N, and therefore each population 
cannot exceed N, despite its fluctuations. However, sometimes it is useful to relax this constraint and allow each 
population to fluctuate within the interval [0, oo). In this case, a parameter may or may not exist that allows one to 
recover the deterministic limit. However, if we assume that [n] are smooth functions of n, and when n —>■ n -f r^, 
they vary little with respect to a characteristic scale of the model (e.g., the total average population), and then, we 


can treat n as continuous variables and expand ti[n — r^]P(n — r^, t) around n. Starting from the ME (30) and using 
X to indicate the continuous populations, we obtain 


ap(x,t) ^ ^ [A,(x)p(x,t)] + J , 


dt 


2 dxidxj 

i,]=i ■> 


(36) 


where time has not been rescaled and Ai(x) and Bij{x) are as in Eqs. ( |33[ ), (34). If we are within a neutral framework 
with one-step jumps, then S' = 1, L = 2 and ri = I, r 2 = —1. Setting ti[x] = h(x) (birth rate) and t 2 [x] = d{x) (death 


rate), from Eq. (36) one gets 


dP{x, = A ^d{x) - b{x)]P{x, 0 + ^ [d{x) + b{x)]P{x, t) . 


dt 


2 dx ‘ 


(37) 


Starting from the discrete formulation of the birth and death rates that we have seen previously, one can derive the 
expressions of the corresponding continuous rates. Eollowing our earlier discussion, we can write the following general 
expansion 


X 



X 


d{x) 

X 



X 


(38a) 

(38b) 


where we have dropped higher order terms in 1/x. Clearly, Eq. (38) can be applied when species have sufficiently 
large populations. The constants &o and do produce a density dependence effect that causes a rare species advantage 
(disadvantage) when bo > do {bo < do). This effect has its roots in ecological mechanisms, such as the Janzen-Connell 
[4ll|29| or Chesson-Warner effects [36| . However, the presence of net immigration or speciation in a local community 
m can also produce such density dependence, which is captured here by a mean field approach. Usually, the 
skewness of the RSA of various tropical forests and coral reefs indicates a rare species advantage |147j , and thus, we 
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will use bo > do- In addition, to simplify the analytical treatment and for parsimony, we choose 69 = —do- Substituting 


(39) 


the rates in (381 into Eq. (371 and setting bo = —do > 0, we obtain 

where Z? = (di + bi)l2, /r = 1/r = di — 61 > 0 and b = 2bo > 0. P{x,t) = P{x,t\xo,0) is the probability density 
function of finding x individuals in the community at time t, given that at time t = 0 there were xq. Therefore 
P{x,t)dx is the fraction of species with a population between n and n + An at time t. Eq. (39) defines our 
model for species-rich communities and it governs the time evolution of a community population under the neutral 
approximation. The deterministic term drives the population to the stationary mean population per species b/^ and 
therefore, within our model, the population and the number of species can also fluctuate at stationarity with the ratio 
of the population to the number of species being fixed on average. Thus, this model is more flexible than the original 
model by Hubbell m in which ‘zero-sum dynamics’ that fixes the total number of individuals in a community was 
assumed. 

The link established between the FP equation and the ME provides a useful interpretation of the coefficients: 
/i = di — 61 is the imbalance between the per capita death and birth rates that inexorably drives the ecosystem to 
extinction in the absence of immigration or speciation. In this model r is the characteristic time scale associated 
with perturbations to the steady-state. When Dr ^ 1 (as for the tropical forests we have analyzed) we have 
di = {2D + t~^)/2 ~ D, and so D can be thought of as an individual death rate. Finally, b plays the role of an 
effective immigration (or speciation) rate that prevents the community from becoming extinct. 

As expected, different choices of bo and do lead to very good fits of the RSA of various tropical forests. In 
particular, it is possible to see that the RSA fits are readily improved for large x when bo and do are arbitrary 
parameters. However, if one introduces a fourth free parameter, the analytical treatment of the dynamics is much 
more involved. The Fokker-Planck equation with just three parameters is an ideal starting point to understand the 
dynamics governing species-rich ecosystems in a simplified fashion. Indeed, the agreement between empirical data 
and the macro-ecological properties predicted by the model is consistent with the simplifying assumption that tropical 
forests are close to their steady state. In addition, because the model is not only neutral, but also non-interacting, 
one may speculate that surviving species are those that were able to reach a steady state of coexistence by minimizing 
interspecies interactions. 


Further insights into the nature of stochasticity can be achieved by writing down Eq. (36) in the equivalent Langevin 


form, which is an equation for the state variables themselves. Within the Ito prescription [83], Eq. (36) is equivalent 
to the following Langevin equation 

s 

Xi{t) = A,(x) -I- ^ (x)^j (t) , (40) 

where ^j{t) is Gaussian white noise with a zero mean and correlation {^i{t)^j{t')) = 6ijS{t — t'), and the matrix &(x) 


is defined by the equation i?(x) = 6 (x) 6 ^(x), where the entries of i?(x) are given in Eq. (34). This shows that b(x) 
may be thought of as the “square root” of B{x). However, because any orthogonal transformation of the noise ^j{t) 
leaves the mean and the correlation unchanged, the matrix 6 (x) is not uniquely determined. 


For the present model, the Fokker-Planck Eq. (391 is equivalent to 

x{t) = b- -f ^Dx{t)^{t) 

T 


(41) 
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where we have instead assumed that the time correlation is = 2(5(t — t'). Eq. (41) shows how demographic 


stochasticity accounts for random ecological drift. In fact, given the relatively large number of individuals of any 
species, one expects that the detailed nature of the stochastic noise is not important, such that fluctuations are 
proportional to x{t) due to the central limit theorem |83j . Thus, the multiplicative noise in the Langevin equation 
is not introduced ad hoc and can be justified on the basis of more general considerations. 

B. Boundary conditions 


The stochastic dynamics described by Eq. (41) governs the evolution of the population when it is strictly positive. 


Assuming that & > 0 and starting at some positive value, the process cannot reach negative values because the drift 
brings it to a positive value bdt in the following time interval dt. Therefore, the process is always non-negative with 
the origin acting as a singular boundary, which must be specified on the basis of ecological considerations, i.e. one 
should define what happens when the random ecological drift occasionally leads to a vanishing population. The most 
frequent ecological situations are given by the following boundaries 

• We use reflecting boundary conditions at a; = 0 to describe ecological communities in which a vanishing pop¬ 
ulation is immediately replaced by a new individual (belonging to either a new or old species). For instance, 
this may happen when an ecosystem is coupled to a met a-community that extends across large spatial scales. 
With these boundary conditions, one can describe ecosystems in which biodiversity is continuously sustained 
by the net immigration of new individuals of new/old species. This prevents an ecosystem from becoming ex¬ 


tinct and can finally achieve a non-trivial steady-state. In the case of the model defined by Eq.(39), reflecting 


boundaries are obtained by setting the flux of the probability distribution as zero at a; = 0, although there are 


some equations for which such boundaries cannot be arbitrarily set, such as Eq.(39) m- 


• A second possibility arises when a community can lose individuals without any replacement or a net emigration 
flux of individuals from the ecosystem exists. One can then describe the system by introducing absorbing 
boundary conditions at a; = 0. These constraints force ecosystems to march inexorably to extinction and 
the final steady-state corresponds to the complete extinction of species, i.e. Pix,t) = d{x). Usually, 

absorbing boundaries are obtained by setting the probability distribution at a; = 0 equal to zero, yet some 


equations do have norm decreasing solutions, such as Eq. (39), m which do not vanish at a; = 0. 


It is possible to define other kinds of boundaries, according to different ecological behaviors when the population 
consists of just a few individuals. However, the reflecting and absorbing boundaries capture the most interesting 
situations in ecology. 


As we have alluded to above, the reflecting and absorbing solutions of the Fokker-Planck equation (39) have rather 


unusual properties. This is essentially due to the fact that the diffusion term, Dx, vanishes at a: = 0. Obviously, 
it is always possible to rewrite a vanishing diffusion term at zero into a non-vanishing one with a suitable change 
of variables, yet the same change of variables leads to a singular drift term at the boundary. Thus, equations with 
a vanishing diffusion term or with a singularity of the drift term at the boundary in general share similar unusual 


properties. In fact, there exist regions in the parameter space of Eq. (39) in which one cannot arbitrarily fix a 


boundary at zero, and the solution is completely defined once the initial condition is prescribed. This is because the 
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probability of accessing the origin depends on the values of the parameters. The problem can be studied by introducing 
the probability to reach the origin for the first time. It is possible to show [96] that this first-passage probability to the 
origin is zero when b/D > 1 and therefore, in this case one cannot obtain any solution with absorbing boundaries and 
the solution can only be reflecting at a; = 0. However, a (unique) norm decreasing solution exists when 0 < b/D < 1, 
which corresponds to the regular absorbing solution but that does not vanish at x = 0. Also, in the same region of 
parameters, a reflecting solution exists with an integrable singularity at x = 0. All ecosystems that we have studied 
so far are well described by b and D within this parameter region. 


C. Stationary and time dependent solutions 


An ecosystem described by Eq. (391 can reach a non trivial steady state when setting reflecting boundary conditions 
at X = 0. In this case, the normalized stationary solution is the gamma distribution [83j . i.e. 


Po{x) = 


_ 2__1 _ 

-x^ e ^ 


(42) 


Tib/D) 

where r(x) is the gamma function [90|. Note the distribution exists only when 6 > 0, i.e. when the ecosystem has 
a net immigration rate of species. Po{x) is the pdf of the relative species abundance at stationarity and it gives the 
probability that a given arbitrary species present in the ecosystem has x individuals. Fitting the empirical RSAs 
allows one to estimate two combinations of the parameters, i.e. b/D and Dt, with r being the characteristic temporal 
scale to approach stationarity. Interestingly, for 0 < b/D 1 one obtains the continuum approximation of the 
celebrated Fisher logseries |52|, the RSA distribution of a meta-community. 


One can analytically solve Eq. (39) at any time with arbitrary initial conditions, so that the evolution of the 


population can be traced even far from stationary conditions. The time dependent solution with reflecting boundary 
conditions {b > 0) and initial population xq is 


Prix,t\xo,0) = ( ^ 




[ipyxoxe 


b 

2D 


1 - e-*P 


X exp 


;^(x-fxo)e 

T , 

■^^XQxe-^P 

1 - e-*P 

D ^ 

1 - 


(43) 


where I^iz) is the modified Bessel function of the first kind [90]. Because Iv{z) — jViv -|- 1) when z —>■ 0+, for 

large time intervals Eq. (43) converges to Po{x) in Eq. (42). By setting x = Dy^/A {y > 0), when b/D = 1/2 the 
process can be mapped into the simpler Ornstein-Uhlenbeck process |83j with a boundary (reflecting or absorbing) at 
X = 0, and when b/D = 1 it can be mapped into the (reflecting) Rayleigh process [54] . 


From the Fokker-Planck Eq. (39), it is easy to derive the evolution of the mean population per species. It is simply 


(x(t)) =bT + (xo — bT)e (44) 

where (x(0)) = Xq is the initial population per species in the community. The mean number of individuals per species 
converges to br at stationarity, with standard deviation ry/bD. The auto-covariance, K{t) = (x(t)x(O)) — (x(t)) (x(0)). 
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under stationary conditions is given by 

poo 

K{t) = / [xxq - {x{t)) (x(0)} ] Pr(x, t\xo, 0)Po(a;o)da;da;o 
Jo 


( 45 ) 


where Pr{x,t\xo,0) is the reflecting solution in Eq. (43) and Po(a;o) is the stationary solution in Eq. (42). However, K{t) 


can be more rapidly obtained from Eq. (41). In fact, since in the Ito prescription ^{t) and x{t) are uncorrelated, one 


gets (^\/x{t)^{t)J = {^■\/x{t)j = 0 because we have assumed that (C(t)) = 0. Therefore the evolution equation for 

{x{t)x{G)) becomes dt {x{t)x{G)) = b{x{0)) — {x{t)x{0)) jx, whose solution is {x{t)x{Q)) = (6t)^ + Hence, 

under stationary conditions, «(!) = bDT‘^e~*^'^ and therefore, the time correlation function is . This shows that r 
is the characteristic relaxation time of the process, the temporal scale upon which the system at stationarity recovers 
from a small perturbation. 


D. The species turnover distribution 


According to the NT, the turnover of ecological communities reflects their continuous reassembly through immigra¬ 
tion/emigration and local extinction/speciation. Species’ histories overlap by chance due to stochasticity, yet their 
lifetimes are finite and distributed according to the underlying governing process. Non-trivial stationary communities 
are reached because old species are continuously replaced by new species, bringing about a turnover of species that 
can be studied and modeled within our framework. 

To measure species turnover one usually considers the population of a given species at different times, then studying 
the temporal evolution of their ratios. Eor an ecosystem close to stationarity, one can look at the distribution P(A, t), 
i.e. the probability that at time t the ratio x{t)/x{0) is equal to A, where x{t) and a:(0) are the population of a given 
species at time < > 0 and t = 0, respectively. Thus, the species turnover distribution (STD) V{X,t) by definition is 


ViXjt) = {6{X - x/xo)) 


= / dxo / dx P{x,t\xo,0)Po{xo)S{X- x/xo). 
Jo Jo 


(46) 


Here P(x,t|a;o,0) is either the reflecting or absorbing solution defined in Eqs. (43) and (49), and Po{x) is given in 


Eq. (42), and A > 0 and t > 0. 


In Eq. (46), one can use the time-dependent reflecting or absorbing solution according to different ecological 
dynamics. We should use the reflecting solution when not concerned with the extinction of the species present at the 
initial time point and especially, when accounting for any new species introduced through immigration/speciation. 
The expression for the reflecting STD can be found in Ref. [T^ . 

One can show that it has the following power-law asymptotic behaviour for a fixed t 


PiX,t) « for A»1 

V{X,t) ^ k 2 {t)X^-^ for A<1 


(47) 


where the functions fci and k 2 are independent of A. Customarily, rather than the random variable A, ecologists study 
r = log(A) that is distributed according to g{r,t) = e^V{e^,t) and that can be compared to empirical data Eigj^ We 
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FIG. 4. [From |12|]. STD for the interval 1990-95 in the BCI forest. The main panel shows the results for individuals of 


more than 10 cm d.b.h., and the inset the results for individuals of more than 1 cm d.b.h. (Center for Tropical Forest Science 


website). The black line represents the analytical solution given by Eq. (46l. 


obtain an estimate of r, the characteristic time scale for the BCI forest, which is around 3500 ± 1000 yrs (for trees 
with > 10 cm of stem diameter at breast height (dbh)) and 2900 ± 1100 yrs (for trees with > 1 cm dbh), where the 
broad uncertainty is due to the fact that the data are sampled over relatively short time intervals. 

These fits not only provide direct information about the time scale of evolution but also, they underline the 
importance of rare species in the STD. b/D is closely tied to the distribution of rare species and in fact, Poix) ~ 
for X Dt. Alternatively, the dependence of the STD on b/D and r suggests that at any fixed time t > 0, rare 
species are responsible for the shape of the STD. 


E. Persistence or lifetime distributions 

A theoretical framework to study and analyze persistence or extinctions of species in ecosystems allows one to 
understand the link between environmental changes (like habitat destruction or climate change p51H51IMlll37lll38j ) 
and the increasing number of threatened species. 

The persistence or lifetime r of a species is defined as the time interval between its emergence and its local extinction 
(see Fig. within a given geographic region (see [5511113| ). In statistical physics, this is known as the distribution 
of first passage times to zero of the stochastic processes describing the species abundance dynamics. According to the 
neutral theory of biodiversity, species can span very different lifetime intervals and thus, at a local scale, persistence 
times are largely controlled by ecological processes like random drift, dispersal and immigration. 


1. Discrete Population Dynamics 


The simplest baseline model to study the persistence time distribution is a random walk in the species abundance 
fi, i.e. ME (101 with b^ = dn = c and absorbing boundary condition in n = 0. According to this scheme, local 


extinction is equivalent to a random walker’s first passage to zero and thus, the resulting persistence-time distribution 
has a power-law decay with exponent 3/2 [5511113j . 


A further step in modeling life-time distributions can be made by taking into account birth, death and speciation 
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FIG. 5. [From [135] ]. Schematic representation of persistence time (or lifetime) of a species r and survival times Ts, defined 
as the time to local extinction of a species randomly sampled among the observed assemblages at a certain time T. 


[SI IlSl 11061 1148j through a mean field scheme of the voter model with speciation |32l [45] (see description in section 
i.e. ME with birth and death rates given by eqs. ( |20| and ( pTj ) in the large Jm limit, b{n) = (1 — i^)n/ Jm and 
d{n) = tl/Jm for n > 0. The corresponding persistence or life-time distribution is given by Prit) = — , where 

P{0,t) is the probability that a species has a zero population at time t. The asymptotic behavior of the resulting 
persistence-time distribution (i.e. Prit)) exhibits a power-law scaling modified by an exponential cut-off |113j : 


Prit) oc t 


— 2^ — L't 


(48) 


for t greater than some lower cut-off value. Here, t is measured in units of Jm and ly is now a speciation rate rather than 
a probability (Eq. (48) is derived in appendix [^. We note that, in this context, Prit) has a characteristic timescale 
Ijv for local extinctions determined by the speciation/migration rate. The general case when none of the coefficients 


is zero and are given by Eqs. (38) can also be solved. In this case, Pr displays a crossover from the t to the t ^ 


behavior at a certain characteristic time and finally, an exponential decay beyond yet another characteristic time scale 
m- It has been shown numerically |19[II34) that for the spatial voter model in dimensions d = 1, 2,3, the exponent 
of the power-law in Eq. (48) changes to t~“, with a <2 depending on the topological structure underlying the voter 


model. In particular, a = 3/2 in d = I and a = 2 for any d > 3, whereas in d = 2 one gets Prit) oc t Inf as 


shown in |115j and references therein. 

We have seen in section jll] that the RSA pattern does not depend on the biological details of the ecosystem under 
analysis. Thus, one may wonder if the persistence time distribution is also a ’’universal” macro-ecological pattern. 
Indeed, it has been shown [I^ll34j that the power-law with an exponential cut-off shape predicted for the persistence 
time distribution by the NT (see appendixis common to very different types of ecosystems (see Figure]^ 

Another interesting and related quantity is the survival distribution Pr^ defined as the probability that a species 
randomly sampled from the community at stationarity is still present in the community after a time t. This quantity 
depends on the initial conditions as Pr^it) = JPr(t\no)Po(no)dno, where Pr{t\no) is the lifetime distribution for a 
species that initially has a population ng. Assuming that the stationary distribution of population abundances is the 
Fisher log series, then Pr^it) ~ t~^ when t ^ t*, whereas Pr^it) ^ when t ^ t* |113[ 1134] . It can be shown 

that this asymptotic behavior of the survivor distribution is valid regardless of the functional shape of the birth and 
death rate b{n) and d{n) involved in the ME driving the evolution of P{n,t) |I35j . 
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FIG. 6 . [Modified from | 134| 1. Comparison between persistence empirical distributions for (a) North American Breeding birds, 
(b) Kansas grasslands, (c) New Jersey BSS forest, (d) an estuarine fish community and the corresponding theoretical species 
persistence times pdfs. The circles and solid lines show the observational distributions and fit, respectively. The finiteness of 
the time window ATu, imposes a cut-off in the maximum observable persistence time and thus, only lifetimes where r < AT-u, 
have been considered and the theoretical predictions have been adjusted appropriately (appendix [b|. 


2. Continuum limit 


In the continuum limit, a crucial distribution for the analysis of species’ extinction is the time dependent solution 


of Eq. (39) with absorbing boundaries at a; = 0. We refer to m for its complete derivation. This probability 
distribution only exists when b < D and is given by 


Pa{x, t\xo, 0) = exp 


X 

Xo 


-e*/’ 




1 - e-*/^ 

fAjr _ I 


(49) 


Note that (49) is finite at a: = 0 but lim 3 ._>o+ Pa(,x,t\xp,Q) ^ 0. 


The lifetime distribution can be calculated analytically using more sophisticated methods through Eq. (49), i.e. 


PT{t) = — /p Pa{x^t\l^Q)dx. It can be shown that the lifetime distribution calculated in this way displays the 
same asymptotic behavior as (48), although the functional form is different |12j . In the continuum limit, one can also 

-1 x-^/D _ I 


obtain the mean extinction time[ 


[ X ‘ 


— X 


-da; = —t (7 -|- ' 0(1 — /?)) 


(50) 


which depends only on 0 < b/D < 1 and r. Here, 7 = 0.577... is the Euler constant and tf}{z) = r'( 2 )/r( 2 :) is the 


logarithmic derivative of the Gamma function |90j . Eq. (42) can be used to fit the RSA of various tropical forests 


yielding (t) /t = 1.94, 1.67, 0.67, 0.95 and 1.38 for Yasuni, Lambir, Sinharaja, Korup and Pasoh, respectively (see 
m)- These time scales are in accord with the estimates of extinction times presented elsewhere [114] and it is quite 
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interesting that (t) jr depends on hjD only, which can be calculated from the steady-state RSA without the need for 
dynamic data. The values of hjD obtained from various tropical forests [12] suggest that (t) ~ r, athough this is not 
built into the model. In general, if r ^ (t), extinction would be much faster than recovery and the ecosystem will 
not reach a steady state. However, if t ^ (t) the ecosystem would recover from external disturbances very rapidly 
with respect to the extinction time and therefore, it would be very robust. This would leave little room for the action 
of evolution. Therefore, the fact that (t) ~ r suggests that ecosystems at stationarity might be marginally stable — 
not so stable that they are frozen in time and not so fragile that they are prone to extinction. From estimates of the 
model parameters h/D, t and thus predictions of (t), many biological and ecological features of the ecosystem may 
be understood. 


IV. NEUTRAL SPATIAL MODELS AND ENVIRONMENTAL FRAGMENTATION 


Space is an essential element to understand the organization of an ecosystem and most empirical observations are 
spatial. Understanding the features and dynamics of an ecosystem cannot, therefore, be disentangled from the spatial 
aspects. Spatial models are basically analytically intractable and the main difficulty is due to the fact that spatial 
models are not equilibrium models ESj. 

Spatial effects can be incorporated into the theoretical framework, with the ME of § remaining formally the same 
by considering the index i in as a composite index i = (a, r), where a identihes the species and r identifies its 
spatial location. The set of all spatial locations r will be denoted by A. Now the transition rates should take into 
account dispersal from nearby locations. To simplify, we will use the notation nj(r), where i indicates one of the S 
species and r a spatial location. 

At present, no coherent spatial neutral theory exists but rather, there is a collection of models and techniques that 
can explain some spatial patterns. In this section, we review some of those approaches. 

The relationship between the number of species and the area sampled is probably one of the oldest quantities studied 
in ecology [153j . Schoener |131j referred to it as “one of community ecology’s few genuine laws”. The Species-Area 
relationship (SAR) is defined as the average number of species (S'(A)) sampled in an area A. 
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FIG. 7. Tri-phasic shape of the Species-Area relationship. The left panel shows the three behaviors on different scales. At a 
local scale the relationship is linear, becoming a power-law relationship at the regional scale and returning to linear at very 
large intercontinetal scales. The right panel shows the same trend from empirical data insi. 
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Arrhenius [S] postulated a power-law relationship {S{A)) = cA^. Empirical curves shows an inverted A-shape |119j 
(see Fig. [^, with a linear behavior at small and large areas, and a power-law with an exponent z at intermediate 
scales. This behavior seems to be pervasive and has been reported for distinct ecosystems. 

Despite some notable exceptions [58], the value of the exponent z has attracted most attention in studies of SAR. 
This value of z is far from universal, ranging from 0.1 to 0.5 [94] and showing dependence on latitude [131] . body-mass, 
taxa and general environmental conditions [9511117] . The exponent z is interpreted as a measure of biodiversity. 

Several models tried to reproduce the empirical behavior and a simple but useful assumption involves considering 
different species as independent realizations of the same process [39] . It is important to note that this assumption is 
stronger than neutrality, because neutrality does not imply independence. Under these assumptions the SAR is given 
by Eq. 0. 

The Endemic-Area relationship (EAR) is defined as the number of species that are completely contained (i.e. 
endemic) in a given area (see BOX 1). It quantifies the number of species that become extinct when a portion of 
landscape is destroyed. It is not generally simply related to the SAR m- If the species are considered as independent 
realizations of a unique process, then 


(E(A)) = StotPoiA^) ■ 


(51) 


where Stot is the total number of species in the system, Po{A^) is the probability that a species is not present in A'^, 
the complement of A, and that it is therefore completely contained in A. 

The /?-diversity is a spatially explicit measure of biodiversity. A simple and useful measure of /I-diversity is the 
similarity index, i.e. the fraction of common species shared between different locations. It can also be defined as the 
probability F{r) that two individuals at a distance r are conspecific [35]. This quantity may be related to the two 
point correlation function Gij(r). Under the assumption of translational invariance, this latter is defined as 


Gy(r) := (ni(x)nj(y + r)) = ^ ^'^'^n,{x)nj{y)6,jS{\\x - y\\ - r) 


SlV 


(52) 


X y 


where ni(x) is the number of individuals of the species i in the location x, | |a: — y| | is the distance between x and y, V 
is the number of site locations and ^(||a; — y\\ — r) selects only pairs at a distance r. F{r) is the probability that two 
individuals at a distance r are conspecific, i.e. it is the ratio between the number of pairs of individuals belonging to 
the same species at a distance r and the total number of pairs of individuals at a distance r. Therefore we obtain 

'EiExEy ni{x)n^iy)6{\\x - y\ \ - r) 


F{r) = 


that can be rewritten as 


F{r) = 


J2x J2y niix)nj{y)S{\\x - y\\ - r) 
'Ei {n^{x)ni{x + r)) _ EiG^^ir) 


EiEj{n^{x)nJ{x + r)) E^EjG^j{r) ' 


(53) 


(54) 


If the spatial positions of different species are independent, then Gij(r) = (n(a;))^ := and by defining Gii(r) = 
G 2 (r), the /3-diversity reads 

G2(r) 


F{r) = 


(55) 


(5-I)p2+G2(r) ■ 

Most theoretical work consists of attempts to relate these ecological quantities with other spatial and non-spatial 
observable factors. For instance, a typical problem is to calculate the SAR knowing the /3-diversity and the RSA over 
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a global scale. One of the future challenges will be to relate and predict spatial patterns on different scales (upscaling 
and downscaling), having only local information on one or more patterns. 


A. Phenomenological models 


Phenomenological models do not assume any microscopic dynamics but they are rather based on a given phe¬ 
nomenological distribution of individuals in space. 

The simplest assumption is to consider individuals at random positions in space [SSI SO]. 

This null model, usually known as random placement, can be used to obtain predictions for the SAR and EAR 
having the RSA or the SAD as an input. Even though this assumption is not realistic, random placement turns 
out to be very useful to capture the relevant aspects of the relationship between RSA and SAR, and it also gives 
reasonable predictions that can be benchmarked against empirical data. In addition, this assumption also allows 
direct connections between SAR and EAR to be formulated iznj. 

Consider S'(Ao) species in a region of total area Aq. Species i has an abundance Ui and the N = individuals 

are uniformly distributed at random in the area Aq. If we observe a sub-region of area A, the probability of observing 
a particular individual is simply A/Ag, while the probability of not observing it is 1 — A/Ag. The probability of not 
observing species i will then be (1 — A/Ao)”b given the fact that the positions of individuals are independent. We 
can then obtain the average number of species observed in an area A as 

(^(A))=5(Ao)-^h--j , (56) 


where S'(Ag) = Stot is the total number of species in the system. 


The simple framework of the random placement model also allows the EAR to be calculated. Using (511 we obtain 


(£(A)) = E(:v 


A 


(57) 


Despite the simplicity of the approximation, the EAR evaluated using random placement captures the quantitative 
behaviour of several observed ecosystems m- 

Under random placement assumptions, one can obtain the EAR from the SAR and vice versa. To calculate the EAR 


in (511, we calculated the number of species with zero individuals in the area complementary to that of interest. This 
number is equal to the difference between the number of species in the whole area and the number of species in the 
complementary area. The complementary area has a non-trivial shape and under general assumptions, this quantity 
is not easy to calculate. Under the random placement assumption, the number of species in the complementary area 
is the SAR of the complementary area. The EAR is therefore. 


(U(A)) = 5(Ao) - S{Ao - A) . (58) 

In [70] . a careful analysis of the reliability of this extrapolation was presented to predict the empirical EAR and it 
was shown that the random placement approximation describes the empirical data well. 
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One can obtain a closed form expression for the EAR and SAR by starting with a RSA distribution. In the case of 
a Fisher log-series (see eq. Q) , the SAR reads 


OO TL / /I \ ^ 

(5(A)) = 5(Ao)-^0-(l-—) = 

which follows from the observation that S'(Ao) is equal to —01og(l 
performed for the EAR, obtaining 


6»log(l + 


r A 
1 - r Ao 


— r) of eq. (18). 


) , (59) 

The same calculation can be 


{E{A)) = -6»log(l - r^) . (60) 

In real ecosystems, individuals are not distributed uniformly in space but rather, due to dispersal limitation, 
individuals of the same species tend to be clustered. This is confirmed by empirical /3-diversity plots |32j . The 
phenomenological approach of random placement can be generalized to a non-uniform distribution of individuals in 
space, taking into account the empirically observed spatial clustering of conspecific individuals [6l]. Individuals are 
distributed in space via a Poisson Cluster Process (PCP). In a PCP, the centers of the cluster of points are uniformly 
distributed in space. A random number of points are distributed around each center according to a given spatial 
kernel. The process depends on two distributions: the number of points in each cluster and the spatial kernel. It 
is possible to show [ST] that these two distributions may be related to the RSA and the /3-diversity. An analytical 
formula for the SAR and the EAR can therefore be obtained given the RSA and the /3-diversity. This approach 
reproduces the ^'-inverted shape observed in empirical systems, showing that this shape can be explained simply in 
terms of spatial correlations of conspecific individuals. 


B. Spatial Stochastic Processes 

There are two possible ways to include space in a neutral stochastic model, either implicitly or explicitly. 


1. Spatially implicit model 


Spatially implicit models are based on the observation that one can relate the sample area A to the total number 
of individuals J m, given that J = pA, where p is the density. One can therefore obtain species-area curves in 
non-spatial models, looking at the scaling of the number of species S with the number of individuals J. Spatial 
implicit models can thus be thought of as mean-field “well-mixed” models, where one can neglect dispersal limitation. 

In a meta-community the number of species ((/)(n)) of a population n is a Fisher log-series (see sec. [ll| 

j,n 

[(j){n))=6—, (61) 

' ' n 


where r is the ratio of the birth to death rate. The average number of individuals in the meta-community is Jm = 
9r/ {1 — r) and the expected number of species can be easily computed 


n—1 


(S'(Jm)) = -6'log(l-f r) = 6»log(l-f ^) . 

77 , H 


Jm ^ 


e 


(62) 


This result corresponds exactly to that found using the random placement in (60). At small sample sizes, the number 


of species Jm 9 is equal to the number of individuals {S{Jm)) = Jm- In other words when small areas are sampled 
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the individuals belong to different species and the number of species grows along with the number of individuals. 
With larger sample sizes, the number of species Jm S> d, grows logarithmically with the number of individuals. This 
approach allows one to calculate the SAR directly from the RSA distribution and it is clearly applicable to any RSA 
distribution. 

Real ecosystems are of course spatially explicit, but one might wonder how spatial implicit models or, more generally, 
models that do not consider space explicitly, are predictive and how their parameters are related to spatially-explicit 
ones. A way to assess this is to measure the efficacy of non-spatial models in predicting the behavior of spatially 
explicit models [SO]. As expected, non-spatial models have a good predictive power when the dispersal lengths are 
sufficiently large and they are particularly good in predicting non-spatial patterns such as the RSA. 


2. Spatially explicit model 


Spatially explicit models are typically defined as birth-death-diffusion processes. A model is fully specified given 
a ME and can be obtained in several ways, i.e. it is possible to write several different MEs that include space in 
a neutral model. The ME is not tractable analytically and one has to introduce approximations in order to get 
analytical results. Spatially explicit models are particularly difficult to solve because of the lack of detailed balance 
(see sec. 0 - 

The voter model m was originally introduced to describe opinion formation, whereby voters are located in a 
network and each one has one opinion among q possibilities. In ecological applications, voters become individuals 
and their opinion corresponds to the species they belong to |45j . An individual (voter) is chosen at random and is 
replaced by a copy of one of its neighbors (see Figure [^. This process has an absorbing state, because when a species 
disappears there is no way to introduce it again. In the long run, the system is populated by only one species. In 
order to overcome this problem, one can introduce the possibility of new species entering the system 1155j , a model 
we will refer to as the multispecies voter model with speciation (MVM). 

Consider a lattice of d dimensions (for a typical ecological landscape D = 2), where a is the lattice spacing with 
exactly one individual at each site (a^ is the average volume occupied by one individual), and a total number of N 
individuals. At each time step, one individual chosen at random is removed and is replaced with a copy of one of 
its neighbors with a probability I — or with an individual of a species not already present in the system with a 
probability ly. The case z/ = I is trivial, whereas the case v — Q has an absorbing state, i.e. a state where a single 
species is present. The MVM is clearly a neutral model, the microscopic dynamics is the same for all species. It is 
also a zero-sum process, because the total number of individuals J is conserved. The mean field version of the voter 
model, where the neighbors of each node are all the other nodes, is governed by the ME Q, with the transition rates 


given by (20)-(21|. 


We want to write an equation for F(r) (the probability that two individuals picked randomly at a distance r belongs 
to the same species). Assuming translational invariance, one obtains 


F(r)‘+i = E(r)‘ (l - ^) + ^ (^(r + e.f + F{r - e.f) , 


(63) 


/i=i 


where F{ry is the probability that two individuals separated by r at time t belong to the same species. The solution 
is obtained with the boundary condition F(0) = I. 
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2D connectivty structure networked connectivity structure 



FIG. 8 . Microscopic moves of the Multispecies Voter Model in a 2—d lattice and with a non-regular network. In the 2D 
structure, each site is occupied by one and only one individual, whose color represents the species it belongs to. At each time 
step, one random individual is replaced by a daughter of one of its neighbors with probability 1 — v (black line). The probability 
that a speciation event occurs is v, wherein the individual is replaced by an individual of a new species (red line). In the case 
of a non-regular network (right panel), the number of neighbors depends on the site considered. 


The stationary solution of (63) can be obtained from a Fourier series in the continuum limit, by taking the limit of 
a —)■ 0, TV —>• oo and —>■ 0, and constraining 7 ^ = 2Dvlo?' to be a constant. Thus, one obtains |155j 


F{r) = 


C7 


1-2 


-{jr) 2 K 2 -D {'jr) . 


(64) 


(27r)‘^/^ 

where Kz{r) is the modihed Bessel function. In one dimension the solution is F{r) oc exp(—r/^), where the length 
^ = 7 “^ is equal to: 


^ = log ■ 


1 - V 


1 - - v)) 

The correlation length depends only on the speciation rate v. It tends to 00 when u 


(65) 

0 , i.e. when speciation 


disappears and individuals of just one species will occupy the entire landscape. The ,5-diversity obtained in (64) was 
shown to be in good agreement with empirical data from different forest censa [32] . The prediction could be improved 
taking into account the Janzen-Connell effect [155] . 

The other important quantity that can be deduced with the voter model is the SAR. It is particularly difficult to 
obtain analytical results, though there are some hints that the behavior can be obtained via a scaling approach [155] . 
which gives an exponent z independent of the speciation rate v. 


Starting from (63) and its stationary solution, one can show that the correlation length scales as 7 “^ or equivalently 
as 1/v^- It is therefore natural to choose as a scaling variable. A second scaling variable can be identihed by 

noting that — (^^)/(’^) ~ Therefore, a second scaling variable can be identified in Using these 

definitions, one can deduce a scaling form for the SAR 

= , ( 66 ) 

where power-law scaling holds when 1 and the exponent z is independent of v. One can also find a scaling 

form for the RSA. Denoting $(A, n, v) as the fraction of species with n individuals in an area A, one obtains 


^{A,n,v)=n . 


(67) 
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By combining the two and fixing the average population per species, one obtains the scaling relationship 

a{2-b) = ^{1-z) , (68) 

where b is constrained to be lower than 2. The exponents can be determined numerically in different dimensions |155) . 
In particular, a reasonable value oi z = 0.3 in two dimensions is obtained. 

It is possible to obtain more precise results via extensive numerical simulations. In the case of NT, one can take 
advantage of the coalescent approach [871 IT29] . Instead of simulating the stochastic dynamics directly, one can 
reconstruct the genealogy of the individuals in the sample area by regressing in time. The main value of this method 
is that there is no need to wait for any transient state to decay and therefore, this approach is much faster than 
forward dynamics [129) . as well as allowing infinite landscapes to be simulated. 

The coalescent approach has been applied to the MVM with a different dispersal kernel |I25j . Instead of a nearest- 
neighbor diffusion, once an individual is removed it is replaced by the offspring of another individual with a probability 
that depends on its distance from the individual removed. In an infinite landscape, the SAR shows the characteristic 
inverted-^" shape. The model depends only on two parameters (up to a choice of the functional form of the dispersal 
kernel): the speciation rate and the dispersal length The SAR scales as 

(5(A,c,^^)) = ^^(A^^^^), (69) 


where the exponent r is independent of v and very close to 2 (as expected by dimensional analysis). The exponent z 
of the power-law regime can be calculated as the derivative of the curve evaluated at the inflection point of the SAR 
(in log-log scale), i.e. the minimum value of dlog(S')/dlog(A). Additional simulations jlllj . obtained with a wider 
spectrum of speciation rates, suggests a logarithmic relationship 


1 


(70) 


q + m\og{v) 

where q and m are two real parameters, which confirms the original prediction of |45j . This inverse logarithmic trend 
seems very robust and it has also been observed in other spatial models m- Indeed, it was shown that power-law 
dispersal kernels better fit data than other short-ranged kernels whereby the estimation of the speciation rate 

corresponding to a given value of the exponent gives much smaller values. 

An attempt to connect spatial and temporal patterns can be found |19j . where the persistence-time distribution 


was studied in MVM (see Sec. HIE). It was shown that the empirical persistence-time distribution is consistent with 
that predicted by MVM and 




(71) 


where t is the average time of persistence in an area A. The exponent a is universal and depends only on the 
dimension of the system, while the time scale r is a function of the area sampled. Empirically, the time scale r scales 
with the area as A^. In the MVM r = Ijv. Using this fact and that the rate of appearance of a new species is 
A = vN ~ vA, one can relate these quantities with the SAR. Indeed, the SAR is the product of the rate of appearance 
of new species and their average persistence time. 


S{A) = \{t) ~ . 


(72) 
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Here we have assumed that (t) ~ to obtain a scaling relationship that connects the exponent z and the other 
exponents: z = 1 — /3(a — 1). This allows spatial patterns to be connected with temporal ones in a manner consistent 
with the empirical patterns shown in Figure 

Other models have been considered in the literature. In a pioneering paper, which received great interest, a ME 
was proposed that seemed to permit analytical calculations m- Indeed, an analytical expression for the SAR was 
obtained. However, the promised solution was not correct and nor was it an approximation of the actual SAR m- 
The main technical difficulty preventing an analytical solution was that the model was inherently out-of-equilibrium. 
Detailed balance, which is the condition that makes the calculation of stationary probabilities possible, did not 
hold [60] . The time is ripe for explorations of this kind. The challenge is to develop new and powerful techniques in 
non-equilibrium statistical mechanics. 


C. Environmental Fragmentation and Habitat Loss 


Resources are not equally distributed and, even over small scales, their distribution in space is far from uniform. 
This spatial heterogeneity affects the distribution of individuals and species in space, and has clear implications for 
the conservation of ecosystems. The space within which ecosystems are embedded is often fragmented and a species 
may be present in only part of the landscape. Moreover, different patches are not independent but they are rather 
connected via immigration. 

In the absence of speciation, the VM m predicts monodominance in dimension D < 2. Spatial heterogeneity can 
be modeled as quenched disorder [23] and when two species are considered, it is postulated that different locations 
on a lattice prefer one species over the other. At each site i, a binary variable ai resides that takes values of ±1, 
depending on which species is present. Spatial heterogeneity can be modeled as a quenched external field which 
also takes values ±1, and the dynamics are fully specified by the transition rates 


W[cri -CTi] 


I - CTiai 
2z 


jedi 


(73) 


where di is the set of nearest neighbors of i, and z is the size of this set. The quantity e measures the strength of 
the preference. The main result from this model [23] was that this randomness enhances the coexistence of species. 
Indeed, if e is larger than = ^J2j(2 N), species coexist in any dimension in the limit of a large system. One 

can introduce the quantity ij) = (l/Af)X)i'^i if coexistence occurs ((()^) < 1, whereas this quantity tends to one 
otherwise. It was shown that [23], for small e 


P{4>) oc 


exp(-f 

(l-e2)(l-()|2) + 2jy ’ 


(74) 


where a small mutation rate {v <C 2/(2 -|- N)) was introduced to regularize the solution. If e > Cc, then the average 
< 1, i.e. coexistence is stable. 

NT has also been applied to predict the extinction rate of species after habitat loss [551 ■ Habitat loss corresponds 
to a reduction of the area available and therefore, to the total number of individuals. When this area is destroyed, 
the endemic species suddenly disappear and what follows is a delayed series of extinctions due to habitat loss. The 
community was modeled as a neutral assembly of species and a typical time scale of extinctions was obtained, along 
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with its dependence on the number of species and the habitat destroyed. The predictions obtained in this way 
reproduced the available data of avifaunal extinctions well [66]. 


V. BEYOND NEUTRALITY AND OPEN PROBLEMS 

A. Reconciling Neutral and Niche Theory 

The concept of niche is central in classical ecology [351123] • An ecological niche is “the requirements of a species 
for existence in a given environment and its impacts on that environment” |29j and it describes how an organism 
or a population respond to changes in resources, competitors and predators. A possible mathematical realization of 
the concept of niche is the Hutchinsonian niche, which involves a n-dimensional hypervolume, where the axes are 
environmental conditions or resources. A position in this space represents a set of behaviors and traits characterizing 
a species or a group of individuals. In niche theory much relevance is given to the specific traits of species and their 
interdependence. A central concept in niche theory is competitive exclusion, which states that two species cannot 
occupy the same niche, as two identical, yet distinct species, cannot co-exist. 

The mathematical representation of an ecological community that includes niche aspects typically coincides with 
the Lotka-Volterra equations. In this case, the focus is on the properties of the fixed points (or other dynamical 
attractors) of these systems of equations and the typical problem that is analyzed is their stability, in relation to the 
parameters and the species present in the system. 

The main difference between neutral and niche theory therefore depends on which mechanism plays the main role in 
shaping ecosystems [ST]. Neutral theory assumes that random processes, such as dispersal, demographic stochasticity, 
speciation and ecological drift, have a stronger impact on many of the observed patterns than niche differences. Niche 
theory assumes the opposite, that the quantitatively important processes are related to differences in species and their 
interdependence. 

As we might expect in a real ecosystem, both stochasticity and niche differences play a role, and it is natural to try 
to quantify how neutral behavior emerges from a niche model. In many cases, niche-based and neutral models yield 
compatible fits of biodiversity patterns [3311331 flOlllIOSj . and it’s impossible to distinguish between the mechanisms 
by looking at those patterns. As pointed out |2] neutrality emerges when species have the same or very similar fitness. 

Neutrality has often been proposed to emerge under some conditions from models considering niche differences m 
1591165j , and neutrality and niche theory were proposed to be the extremes of a more general model [591 165] . In 
both cases, a Lotka-Volterra equation is introduced to describe community dynamics, while population dynamics is 
modeled by also taking into account demographic stochasticity and immigration. By considering different values of 
parameters, one can move from a scenario where species differences matter a lot and the stable configuration is very 
close to the solution of the deterministic Lotka-Volterra equation, to a scenario where demographic stochasticity is 
more important and the community behaves like a neutral community. A slightly different approach has also been 
considered |53j . analyzing a stochastic version of Lotka-Volterra dynamics and quantifying, when stochasticity is 
varied, the difference between the prediction of the neutral model with the full model. In this case, instead of a 
continuum of strategies, neutral and niche regimes are two macroscopic phases separated by a phase transition. 

A different approach was based on effectively considering niche theory as a model where per capita death/birth 
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rates are species dependent [21]. In this case, one might expect the difference between species abundances to reflect 
the differences between these parameters. More precisely, in a neutral scenario, all the species fluctuate around a 
given abundance value, while when niche characteristics play a role, each species fluctuates around its own distinct 
value of abundance. A third scenario has been proposed [24], wherein the per-capita birth/death rates do not have a 
monotonic effect on species abundance. The symmetry between species, due to the neutrality of the process, can yet 
be spontaneously broken. In this case, the stable states are not symmetric in terms of species abundance. 


B. Emergent neutrality 

Ecologists have highly criticised NT because of its unrealistic assumptions. The patterns that we have studied so 
far can in fact be explained without introducing species differences, and this has led some ecologists to oppose NT 
because it assumes that nature is actually governed by neutral processes, whereas it is not. Clearly, no one believes 
that nature is truly neutral. The patterns of community ecology are actually generated by a cocktail of processes, and 
it is both inappropriate and dangerous to consider processes in isolation from a macro-ecological pattern or empirical 
data set. 

However, it is informative to study whether, how and which ecological processes can drive a community towards a 
state in which demographic stochasticity and immigration play a crucial role in the face of strong species’ differences 
that are dictated by classical competitive exclusion. Such a state should allow similar species to emerge in the niche 
space with the ability to co-exist for sufficient time. Indeed, when this problem was studied, it was shown that species 
can evolve into groups of relatively more similar species that co-exist for very long times m- First, a large number 
of species were placed at random along a hypothetical niche axis, which represents a specific trait, assuming that 
interspecific competition can be calculated through niche overlap. Running a classical Lotka-Volterra competition 
model and studying evolution, groups of multiple species were evident that aggregated around similar values in the 
niche axis and they could co-exist for many generations before the majority of them head towards an inexorable 
extinction. Eventually, only one species survives from each group, producing the expected pattern of single species 
equally spaced in the niche space. 

In other words, the niche similarity of species prevents competitive exclusion from swiftly selecting the best com¬ 
petitor among a group of similar species, allowing their co-existence for very long times even though only the superior 
species will ultimately persist. 

This model may be considered as one of the possible steps towards a reconciliation of niche and neutral theories. 
Species that are initially ecologically non-equivalent, and that therefore behave in a non-neutral fashion, are driven 
by community and evolutionary processes towards states in which the dynamics may well be better approximated by 
neutral models over appropriate spatial and temporal scales. More recently, further support for this approach came 
from showing that the model can produce multimodal RSAs m- Immigration may also be an important component 
in the neutral-like behaviour of communities isniiTi]. Parasitoids competing for a common species [22] have been used 
to show that clusters of species separated by gaps emerge along the niche axis, confirming — using a quite different 
approach — that processes exist that can lead community dynamics to be effectively neutral [I20j . 

The basic Lotka-Volterra model has been extended to investigate the possibility that some processes decrease 
the risk of competitive exclusion so that species lumps are not only transient, but ultimately permanent m- 
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Density-dependent regulation was introduced that stabilizes the co-existence of species within a group. This approach 
unfortunately has a drawback that the mechanism introduces a discontinuity in the competition strength of the species 
m, which means that unmodeled species differences may be responsible for co-existence in the community. However, 
it has been shown that more realistic density dependent terms, or even other mechanisms (e.g., migration |144j l. can 
eliminate this problem, making the approach more robust |146) . 


C. Maximum Entropy Models 

The maximum entropy principle is a useful method to obtain the least biased information from empirical measure¬ 
ments m- This powerful tool, borrowed from statistical mechanics, has a wide range of applications m, including 
ecology [niiiT]. In its ecological application, the Max Ent principle is an inference method [34] used to evaluate 
the effective strength of interactions among species based on either species abundance data |150j or simply the pres¬ 
ence/absence of the species |llj . This methodology provides a way to systematically incorporate the most important 
species interactions into the development of a theory beyond the purely non-interacting case. In addition, the Max 
Ent principle was implemented as a method to predict biodiversity patterns across different spatial scales using only 
the information on local interactions |3l|67|. Here we will discuss the extension of Max Ent to the study of spa¬ 
tial biodiversity patterns. Consider an ecosystem in which S species (belonging to the same trophic level, see BOX 
2) live within a given area H, divided into N adjacent sites of equal area. Let us assume that there are empiri¬ 
cal records of the species contained within each site. From these data one can calculate, for example, the average 
presence of any species in the ecosystem and the co-occurrence of any pair of species in neighboring sites. When 
applying Max Ent, we can consider these mean occurrences and co-occurrences as given constraints. We introduce 
the binary random variable cr“, which records the occurrence of species a = I, ...,5" at each site i = I,...,TV. If 
species a is present in plot i, then uf = I, otherwise cr“ = 0. Therefore the “state” of any species a can be char¬ 
acterized by the random vector (t“. There is empirical evidence that species belonging to the same trophic level 
interact weakly [nKnmsoj and therefore, in a first approximation, one may assume that species occur in a given 
geographical location independently of one another. Due to such independence, the probability of finding the system 
in the configuration cr = (cr^, cr^,.., cr'®) is thus P{cr) = where Pa gives the probability distribution 

of finding a species a in the configuration (t“. In order to build the Max Ent model, we will maximize the Shan¬ 
non’s entropy H = — '' 'J2ctS ~ imposing the average occurrence 

constraint, i.e. (M) = J2a Pa{cr°‘)Ma{cr°‘) with Ma{cr°‘) = J2i average co-occurrence constraint, i.e. 

{E) = S<T° Pa(cr“)Lia(<7'“) with Ea{cr°‘) = J2(i j)^ where {i,j) indicates two nearest neighbour locations. 
Both {M) and (E) are meant to match the corresponding empirical averages, as calculated from the empirical records 
of species contained in the region. Maximization |S] provides an explicit expression for the probability, pa, of finding 
a species a in the configuration (t“. Thus, 


Pa{fr ) = 


„ JcSq (<t“)-|-;i„M„ (cr“) 


(75) 


where = Za{ha, Ja) is the partition function. From (75) one can characterize the spatial biodiversity patterns of 


the ecosystem, i.e. calculate the SAR or the EAR |3]. This type of approach may well be suited to infer biodiversity 
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properties of communities over larger spatial scales by upscaling the model results at local scales. However, only a 
few studies have tackled this problem (see below). 


D. Linking different Macro-Ecological Patterns 

Over the last few decades, ecologists have come to appreciate the importance of spatial patterns and processes, 
and the explicit introduction of space has the potential to revolutionize what we know about natural populations 
and communities |13311139j . It has become apparent that key ecological patterns, such as SAR, RSA and spatial 
patterns of species distributions and turnover, are intimately intertwined and scale dependent (see Fig.l). However, 
despite a blossoming of models to address spatial patterns, only a few practical methods have been proposed to 
link them across different scales. Specifically, spatial approaches [10011116] lack the needed analytical machinery, 
whereas most theoretical approaches are not spatially explicit or sufficiently flexible [551 1150] . Even neutral theory 
was conceived to reflect the idealized behavior of natural systems at equilibrium, rather than to reflect non-pristine 
landscapes produced by environmental change or management m- Therefore, a general methodology is required to 
predict and link these ecological patterns across scales that is sufficiently robust and flexible to allow its application to 
a range of natural or managed systems. One possible way to tackle this challenging problem is through a theoretical 
framework inspired by ideas coming from phenomenological renormalization m- The fundamental assumption at 
the core of this theoretical setting is that the functional form of the RSA remains the same across all spatial scales, 
even though the parameters of the curve are likely to vary. Because of this assumption, the spatial dependence of the 
abundance distribution can be obtained by making the RSA parameters suitable functions of scale. Together with the 
functional shape of the RSA, the other model input is the spatial pair correlation function (PCF), which describes 
the correlation in species’ abundances between pairs of samples as a function of the distance between them [155] . If 
populations were randomly distributed in space, distinct communities would on average share the same fraction of 
species regardless of their spatial separation, and therefore, the PCF would not depend on distance. In contrast, in 
highly aggregated communities correlations in abundance would fall off steeply with increasing distance. The PCF 
not only measures the rate of turnover in species composition but it also reflects the variation of population clustering 
across scales, given that the variance in species abundances at any particular scale can be calculated directly from 
the PCF [TU]. Therefore, the PCF is related to the spatial species abundance distribution. Thus, the PCF can link 
the effects of aggregation, similarity decay, species richness and species abundances across scales. Building on the 
intrinsic relationship among these patterns, the profile of the species area relationships and the species abundance 
distributions could be predicted across scales when a limited number of fine-scale scattered samples is available [10] . 


E. Multi-Trophic Ecosystems: Ecological Networks 

Although we are progressing in understanding the suitability and limits of non-interacting and non-spatial models, 
most neutral models still assume that species interact randomly with each other. However, a network approach 
to modelling ecological systems provides a powerful representation of the interactions among species [HI [Hi 1136]. 
Ecological networks may be viewed as a set of different species (nodes) and connections/links (edges) that represent 
interspecific interactions (e.g., competition, predation, parasitism and mutualism). The architecture of ecological 
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interaction networks has become a bubbling area of research, and it seems to be a critical feature in shaping and 
regulating community dynamics and structure diversity patterns [5]. An important step relevant for multi-trophic 
systems will be to obtain a general framework within which a network of preferences/disfavors modifies the birth 
and death rates of different species and can be superposed on neutral models (like the voter model presented in this 
review). 


VI. CONCLUSIONS 

In this paper, we have attempted to describe some of the theoretical frameworks that can be used to understand 
key issues related to biodiversity and that will serve to address important questions. These frameworks are neces¬ 
sarily elementary and incomplete, yet they have the advantage of being tractable and related to the central issues. 
Unlike standard approaches to more traditional physics, here the Hamiltonian function or the interactions among the 
components are completely unknown and even identifying the state variables may sometimes be a non-trivial task. 
We have a fleeting picture of what an ecosystem is and it is not necessarily in equilibrium. We have little knowledge 
of the myriad degrees of freedom and their interactions. The real challenge is to discern the most essential degrees of 
freedom and to develop a framework to understand and predict the emergent ecological behaviour. 

In the near future, one needs to investigate how important the species traits are and their interactions as well as 
environmental effects. At the moment, we do not have a spatially explicit theory that can predict analytically the 
most important ecological patterns. This would be important, because it would allow us to understand what drives 
biodiversity across spatial and temporal scales. Indeed, it is likely that biological processes are not equally important 
across scales. As evident in particle physics, we might eventually find that ecosystems will need to be described by 
different effective theories according to the range of scales in which we are interested. Moving towards even more 
profound questions, a unified theory would ideally explain what determines the total number of species and individuals 
that live in a given region. Throughout this review, S and N are parameters which have to be provided as input. 
Currently, we do not have a mechanistic explanation for why one region may be more bio-diverse than another, what 
sustains biodiversity and how evolutionary pressures sculpt ecological communities. 

The basic message of this review is that to resolve these challenging problems, ideas and techniques must be recruited 
from different disciplines. We are still at the beginning of this adventure. Moving forward is not only important but 
it is also urgent. The pressures of habitat destruction, pollution and climate change are having highly undesirable 
consequences on the health of ecological communities. To address practical issues related to conservation biology, we 
need models that can be used across scales in order to extrapolate information on biodiversity from accessible regions 
to inaccessible yet important scales. There is plenty of room for ideas that matter, and community ecology can greatly 
benefit from the contribution of other disciplines, including physics. 
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Appendix A: Density dependence, and assembling of local communities. 


In this appendix, we present an alternative method to introduce density dependence and how ecosystems emerge 
by assembling local communities. 

First, one can set bk{n) = b ■ {n + T^) (for the fc-th species) and d{n) = dn, where incorporates both the effects 
of intra-specific interactions, such as those giving rise to density dependence, and the immigration occurring from a 
met a-community. This approach can be applied to two distinct ecosystems: coral reefs and tropical forests. 

The steady-state solution of the ME for Pk{n) yields a negative binomial distribution: 


Pk{n) 


(1 - rk)'^’‘ 

r(Tfc)^r(n + Tfc) 


(Al) 


with a mean (rife) = rfcTf;/(l — Vk), where T is the gamma function. 

The number of species containing n individuals is given by (/>„ = J2k=i ^n,k, where In,k is a random variable that 
is 1 with a probability Pk{n) and 0 with a probability 1 — Pk{n). Thus, the RSA is given by 


((^„) = ^/„,fe = ^Pfc(n)=0^r(n + T), (A2) 

k=l k=l 

where 9 = £'/[(! — r)“^ — l]r[T] is the biodiversity parameter [T^, and we have dropped the k dependence because 
of the symmetric hypothesis. For a small T, the RSA for the communities resembles the Fisher log-series and it does 
not have an interior mode. 

Note that a non-trivial k dependence might arise even under the neutral hypothesis. For example, one can set 
T^, = rhpk, where m is a measure of the immigration rate from the met a-community, in units of the birth rate b and 
Pk is the fraction of individuals in the surrounding meta-community belonging to the fc-th species. As a result, one 
can obtain the following RSA for the community of tropical forests. 


N .a;" /■°°F(y + n) , 

Wn)=9—I ^dy = 9 — f{n,uj), 


(A3) 


nl Jo r(y-f 1) 

where w = A — ln(l — x). This approach provides a virtually indistinguishable fit to the empirical data as (<()„) = 
considered earlier with the advantage of having ecologically meaningful parameters. 

The average number of species observed in the local community is 


Sobs = {S)=S- {<Po) = ^(1 - r)-^ = ^[1 - (1 - r)^] 


(A4) 


fc=i 


If the sample considered has Jl individuals, and thus the community dynamics obeys the zero-sum rule (i.e. it has 
a fixed total population), then the multivariate probability distribution is 

Jl +J2k'^k - ^ 

Jl 


P(n| Jl) = Af PJ Pk{n)S{JL-ni-n2-...-ns) = 


k=l 


n 

k=l 


+ Tfc - I 
nk 


S{JL-ni-n2-...-ns) 

(AS) 


that is a compound multinomial Dirichlet distribution where Af is the normalization constant. 
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Meta-community Composition. Now, let us gradually assemble the met a-community of coral reefs by consider¬ 
ing it as an assemblage of local communities. Let us start by considering the joint RSA of two local communities, A and 
B, with UA and ns individuals, respectively. The probability that the species has n individuals in the meta-community 
formed by A and B is |149) 

j,n 

P{n = nA + nB)= P(n^)P(nB) oc —r(n-I-2T), (A6) 

n\ 

nA+nB=n 

where the actual spatial locations of the local reef communities have been neglected (all local communities are well- 
mixed in the meta-community, i.e. a mean field approximation). The elegant result of the T’s adding to each 
other follows ecologically from their interpretation as immigration rates. For the meta-community, we can introduce 
speciation with a rate ;/ <C 1 as we have seen it has a crucial role when n = 0 (under the assumption of neutrality, 
the species label of the new species is of no consequence), i.e. Tfc = z/ and P{nk) = vr"^ jn P 0{u'^). 

Extending the calculation of the joint RSA distribution to more and more local communities, it can be shown that 
the RSA of the meta-community is characterized by an effective immigration parameter LT, where L is the total 
number of local communities comprising the met a-community, and it becomes log-normal-like if L ^ 1, in agreement 
with the available data |149| . 

A Fisher log-series is observed in two limiting cases - in the met a-community in which there are no immigration 
events and in the very small local community that has a high immigration rate from the meta-community characterized 
by a Fisher log-series RSA. 


Appendix B: Persistence or lifetime distributions 


In this appendix, we present the derivation for the persistence (or lifetime) distribution that leads to the empirical 


choice made in (48). The master equation we want to solve is 
dP{n,t) 


dt 


= b{n — l)P(ji — 1, t) -I- d{n + l)P{n -|- 1, t) — (b{n) + d{n)) P{n, t) , 


(Bl) 


where the birth rate is b{n) = (1 — v)n/JM with b{—l) = 0, and the death rate is d{n) = u/Jm and n > 0. We 
can redefine the time scale t —>■ JMt so that the factor 1 /Jm disappears from the birth and death rates. As the 
initial condition we choose P(n,t = 0) = 6n^no with the initial population as Uq. Our goal is to calculate the survival 
probability defined as 


^(^1^) = ^ Pi.'n, t) = 1 — P{n = 0,t) . 


n>l 


Thus, we introduce the generating function 


G(z,t) = J^P(n,t )2 


(B2) 


(B3) 


n>0 


whose radius of convergence is > 1. Note that G(z = l,t) = 1 and G(z = 0,t) = P{n = 0, t) = 1 — V{t\iy). Using Eq. 

itely 

dG{z,t) 


(Bl I the time evolution of the generating function is derived immediately 

dG{z, t) 


dt 




dz 


(B4) 










46 


with the initial condition G{z,t = 0) = The previous equation is a linear partial differential equation and it can 
be solved by standard methods. One introduces a function Z{t) satisfying the time evolution equation 

dZ{T) 


dr 


= _[(l_^)Z(r)-l](Z(r)-l) 


with a ’’final” condition Z{t = t) = z. Using Eq. (B4|, one is led to 

dG{Z{T),t) 


dr 


= 0 , Vr , 


implying that 


G{z,t) = G{Z{t),t) = G{Z{0),0) = Z(0)’^° . 


The solution of eq. (B5|, with the chosen ’’final” condition gives 


(1 - u)A{z,t) 


1 - (1 -v)z 


and thus. 


G(z,f) = 


1 — A{z, t) 

1 — (1 — v)A{z,t) 


Finally, we get the survival probability 

V{t\v) = 1 - G(z = 0,t) = 1 - (l + ^ - 1)”^) 


- 1 \ “"0 


In the scaling, if we consider the limit of fixed vt as t becomes large, we get the scaling form 

1 


with 


The lifetime distribution is simply given by 


V{t\u) = -F{vt) , 


F{x) = no 


- 1 


Prit) = -dV(t\v)/dt = , 


where the last equality holds in the scaling limit as 

f{x) = noe^ 


— 1 


(B5) 


(B6) 


(B7) 


(B8) 


(B9) 


(BIO) 


(Bll) 


(B12) 


(B13) 


(B14) 


which for X —>■ 0 tends to a non-zero constant. As x —>■ oo, /(x) decays exponentially, leading to the empirical scaling 


form of the kind given by (48). 


Appendix C: The species turnover distribution with absorbing boundary conditions 


In this appendix we present another temporal pattern predicted by the model defined by (39) or (41). For an ecosys¬ 


tem at or near stationarity, the model can provide the exact expression for the STD when using the time-dependent 
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absorbing solution of (391: the so-called the species turnover distribution with absorbing boundary conditions. Unlike 
the reflecting species turnover distribution, this distribution corresponds to a new kind of measure that only accounts 
for the species present at the initial time, and it does not take into account any new species introduced by immi- 
gration/speciation or any old species that reappear after their apparent extinction until t > 0. This amounts to the 
selection of a particular sample and the study of the temporal behavior of those selected individuals. The species 


turnover distribution with absorbing boundary conditions, can be obtained through formula (46) where 


the conditional probability must now be Pa{x^t\xo,0) defined as in (49). The final expression is 




sininb/D) e-‘/^ (e‘/^ - l) 
TT{l-b/D) ^(A+l)2 
/ 3 4Ae“‘/'" 


b/D 


(Cl) 


where b/D < 1. Note that decays exponentially to zero at a rate (1 — /3)/t, regardless of A. 

It is noteworthy that the STD’s with absorbing and reflecting boundaries are indistinguishable whenever t <C r and 
A are not too small (for b/D < 1). Since the BCI data are sampled over relatively short times intervals (at most 10 
years), the distributions are almost the same within the interval 1/2 < A < 3/2. 
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